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ABSTRACT 

As  a  consequence  of  the  ever-shrinking  sizes  of  nanoelectronic  devices,  hitherto 
neglected  quantum  effects,  such  as  tunneling,  are  becoming  important  for  device 
characterization.  The  study  of  electron  reflection  and  transmission  probabilities  at 
potential  barriers  is  one  of  the  important  areas  of  active  research  in  this  field. 

Analytic  solutions  for  the  quantum-mechanical  transmission  coefficient  through  a 
potential  energy  profile  of  arbitrary  shape  do  not  exist.  One  conceivable  method  for 
finding  the  transmission  coefficient  through  such  a  potential  involves  transfer  matrices. 
This  technique  is  numerically  limited,  unfortunately,  and  fails  to  provide  adequate  results 
for  potentials  of  interest  in  the  development  of  practical  nanoelectronic  devices. 
However,  within  its  capabilities,  the  transfer  matrix  method  is  a  useful  reference  to  which 
other  results  may  be  compared.  Another  method,  utilizing  backward  recurrence,  has  been 
proposed  as  a  numerically  stable  alternative  for  calculating  the  transmission  coefficient 
through  such  potentials.  This  second  method  has  yet  to  be  widely  applied. 

This  thesis  investigates  the  capabilities  and  limitations  of  each  method,  with  an 
emphasis  on  their  scope  of  applicability.  Extensive  programming,  in  the  C  language,  has 
been  done  to  examine  the  two  methods.  Output  from  these  programs  has  been  analyzed, 
and  the  backward-recurrence  method  has  been  shown  to  have  wider  applicability,  and  to 
be  faster  and  much  more  numerically  stable. 
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I.  INTRODUCTION 

Numerical  methods  are  required  to  predict  the  characteristics  of  nanoelectronic  devices, 
which  rely  on  quantum  tunneling  for  their  operation.  Continuing  scale  reduction  of  modern 
semiconductor  electronics  will  inevitably  result  in  devices  which  rely  on  the  principles  of  quantum 
physics. 

In  this  thesis,  I  shall  apply  two  different  numerical  methods  to  the  problem  of  calculating 
the  transmission  coefficient.  One,  which  relies  on  transfer  matrices,  is  known  to  have  limited 
applicability  due  to  inherent  numerical  instabilities,  but  provides  a  useful  benchmark  in  cases  for 
which  it  is  accurate.  The  other  is  a  backward-recurrence  algorithm,  thought  to  be  much  more 
numerically  stable,  which  potentially  has  wide  future  applicability.  I  have  written  several  original 
computer  programs,  in  the  C  language,  which  utilize  both  the  transfer  matrix  and  backward- 
recurrence  techniques.  These  programs,  and  analysis  of  their  output,  are  my  contribution  to  this 
research.  To  my  knowledge,  the  backward-recurrence  technique  applied  herein  has  not  been 
widely  used;  it  proves  to  be  powerful  in  solving  for  the  transmission  coefficient  through  difficult 
potentials. 

The  values  for  the  transmission  coefficient  obtained  by  each  method  are  first  compared  to 
a  known  correct  result:  the  transmission  coefficient  through  a  single,  square  potential  energy 
barrier.  In  addition  to  this  test  for  accuracy,  the  transfer  matrix  method  contains  a  built-in  test  for 
numerical  instability,  which  also  is  exploited.  Examples,  which  reveal  the  numerical  limitations  of 
each  method,  are  included  as  well.  Graphical  program  output  is  analyzed  to  critique  the 
capabilities  of  the  two  numerical  methods. 
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II.  QUANTUM  TUNNELING 

In  classical  physics,  the  motion  of  a  particle  is  governed  by  conservation  of  energy:  a 
particle  of  total  energy  E  is  unable  to  move  past  any  point  beyond  which  the  potential  energy  is 
greater  than  E.  Such  a  point  is  said  to  comprise  a  potential  barrier  to  the  particle's  passage.  If 
electrons  in  semiconductor  devices  truly  behaved  as  classical  particles,  then  the  tunneling  diode 
and  zener  diode  would  not  exist,  and  aluminum  household  wiring  could  not  conduct  current 
through  twisted-wire  junctions  covered  by  an  oxide  layer  which,  in  bulk  form,  is  an  excellent 
insulator.  These  items  all  depend  upon  the  concept  of  quantum  tunneling. 

Quantum  theory  allows  for  the  penetration  of  barriers  of  energy  greater  than  the  incident 
particle's  energy.  This  phenomenon  is  a  consequence  of  uncertainty  in  the  position  of  the 
electron,  according  to  the  uncertainty  principle.  (Eisberg  &  Resnick,  1985,  pp.  65-68)  The 
electron's  wave  function,  ^(x,t),  has  a  non-zero  magnitude  at  all  points  in  space,  though  it 
decreases  rapidly  beyond  potential  barriers.  Since  the  probability  of  finding  the  electron  in  a 
specified  region  is  proportional  to  the  square  of  the  magnitude  of  the  wave  function,  there  results 
a  finite  probability  that  the  electron  will  be  found  on  the  other  side  of  the  barrier,  in  a  classically 
forbidden  region. 

For  non-relativistic  systems,  motion  is  governed  under  quantum  theory  by  the  Schrodinger 
wave  equation.  (Eisberg  &  Resnick,  1985,  pp.  177-209)  The  general  three-dimensional  form  of 
this  equation  (the  time-dependent  Schrodinger  wave  equation)  is 

iW         h2 
ih  —  = V2x¥+V(x,y,zy¥.      (1) 

dt         2m 


If  the  potential  of  interest  is  essentially  one-dimensional  in  form,  Equation  1  simplifies  to 

..  ay      h2  d2v   __,  XUJ     „. 

ih = -  +  V(xf¥.       (2) 

dt         2m  dx" 

In  the  above  equations,  >i=h/27t,  where  h  is  Planck's  constant  [6.626xl0"34  J  sec],  /  is  the  square 

root  of  -1,  m  is  the  mass  of  the  particle,  V(x)  is  the  potential  function,  x  is  the  spatial  coordinate, 

and  /  is  time.  In  the  development  of  this  equation,  several  key  hypotheses  must  be  made,  which 

are  worth  reviewing. 

First,  the  wave  equation  must  be  consistent  with  de  Broglie's  postulate  regarding  the 

wave-particle  duality  of  matter,  namely 

'-I-    (3> 

where/?  is  the  momentum  of  a  particle,  and  X  is  the  wavelength  of  the  wave  packet  representation 
of  the  particle.  (Thornton  &  Rex,  1993,  pp.  178-179)  Second,  the  wave  equation  must  be 
consistent  with  the  relation 

E  =  -^  +  V,       (4) 
2m 

the  non-relativistic  expression  for  total  energy  of  a  particle  as  the  sum  of  its  kinetic  and  potential 

energies. 

It  should  be  noted  that,  although  Equation  4  holds  only  for  particles  moving  at  non- 
relativistic  speeds,  experience  has  shown  that  the  Schrodinger  formalism  is  adequate  to  describe 
the  behavior  of  the  systems  studied  in  this  thesis  research.  (Eisberg  &  Resnick,  1985,  pp.  128- 
132) 

The  form  of  the  wave  function,  which  is  the  solution  to  Equation  2,  must  next  be 


determined.  For  a  free  particle,  one  form  of  this  solution  is 

^(xj)  =  cos(kx  -  cot)  +  i  sin(kx  -  cot),       (5) 

where  k  is  the  wave  number  [2%IX\,  cois  the  angular  frequency  [27tf|,  t  is  time,  and/is  the 
frequency.  This,  as  required  above  for  the  case  of  V  constant,  is  a  traveling-wave  solution.  It  is 
also  complex-valued,  which  is  critical  to  the  remainder  of  this  discussion.  The  physical 
significance  of  this  complex-valued  wave  function,  postulated  by  Born  in  1926,  is  that  the 
magnitude  squared  of  the  wave  function  is  the  probability  density  of  the  particle,  or,  stated 
another  way,  that 

^\xjy¥(xj)dx       (6) 

is  the  probability  of  finding  the  particle  between  x  and  x+dx,  at  time  /.    f^fxj)  here  is  the  complex 
conjugate  of  the  wave  function. 

If  the  potential  of  interest  is  not  a  function  of  time,  so  that  the  system's  energy  is 
conserved,  then  the  left-hand  side  of  Equation  2  may  be  associated  with  the  characteristic 

equation  ih =  EY  ,  where  E  is  the  energy  eigenvalue  (a  constant)  of  the  system.  The  spatial 

dt 

and  temporal  dependencies  of  5Pmay  then  be  separated,  so  that  W(x, t)  =  tp(x)e       . 

The  result  of  this  separation  of  variables  is  that  the  time-independent  Schrodinger  wave 

equation  is  given  by 

~^-+V(x)¥(x)  =  EW{x%       (7) 

2m  ax 

where  V(x)  is  the  potential  and  E  is  the  total  energy.    y/(x),  here,  represents  only  the  spatial 
dependence  (on  x,  in  this  case)  of  the  wave  function  T(x,t).  The  lull  derivation  is  given  by  Eisberg 


and  Resnick.  (Eisberg  &  Resnick,  1985,  pp.  15 1-167) 

Equation  7  is  an  ordinary  differential  equation,  not  a  partial  differential  equation,  this  is  a 
major  reduction  in  the  complexity  of  the  problem.  We  shall  now  focus  on  the  solutions 
(eigenfunctions),  y/(x),  of  Equation  7,  which  hold  for  a  region  of  constant  potential  V: 

w(x)  =  Ae,kx  +  Be,kxfor  E>V 
y/{x)  =  Ce "a  +  De^Jor  E<V, 

where  A,B,  C,  and  D  are  constants,  and  k  is  equal  to  ire,  with 


K  =  ^-\ '-.        (9) 

n 

Equation  9  provides  the  correct  propagation  constant  only  when  V  is  constant.  While 

approximating  V(x)  as  a  constant  may  seem  limiting,  it  will  be  shown  that  one  can  obtain  useful 

results  in  spite  of  this  approximation.  Note  the  presence  of  the  increasing  exponential  eK\n 

Equation  8.  This  observation  will  later  prove  critical. 

The  boundary  conditions  for  the  spatial  wave  function  are  that  y/(x),  and  its  first  derivative 
with  respect  to  x,  must  be  finite,  single-valued,  and  continuous.  Eisberg  and  Resnick' s  description 
of  the  reasons  for  these  restrictions  is  excellent.  (Eisberg  &  Resnick,  1985,  pp.  155-157)  We  will 
utilize  these  boundary  conditions  in  what  follows. 

As  an  example,  let  us  apply  these  solutions  to  the  case  of  a  rectangular  potential  barrier. 
Let  this  barrier  potential  be  V(x)=0  for  x<-a  and  x>a,  V(x)=V0  for  -a<x<a.  Classically,  if  the  total 
energy  of  a  particle  is  such  that  E>V0,  then  the  particle  always  passes  the  barrier;  similarly,  if 
E<V0,  then  the  particle  never  passes  the  barrier.  Since  the  spatial  wave  function  if/(x)  has  the 
form  of  Equation  8,  for  the  quantum  mechanical  case,  the  behavior  of  the  particle  is  different. 


First,  i£E>Va,  there  will  be  a  nonzero  probability  of  reflection,  at  both  edges  of  the  barrier. 
Stated  differently,  there  is  less  than  a  100%  chance  that  it  will  traverse  the  barrier.  More 
importantly,  however,  \fEV0,  there  will  be  a  nonzero  probability  of  transmission  through  the 
barrier,  or  tunneling,  also  known  as  barrier  penetration. 

The  form  of  y/(x)  in  the  region  x   -a  is  Ae,kx+Be~,kx,  corresponding  to  a  pair  of  incoming 
and  outgoing  traveling  waves.  The  form  of  the  general  solution  in  the  region  x  a  is  similar, 
namely  Fe,kx+Ge,kx.  Between  -a  and  a,  however,  the  form  is  Ce'KX+De'a.  Boundary  conditions 
require  that,  at  x=-a,  Ae'ka+Be'ka=CeKa+DeKa  so  that  y(x)  is  continuous.  In  order  for  the  first 
derivative  dy/(x)/dx  to  be  continuous,  Aelka-Belka=(iK/k)[CeKa-DeKa].  It  is  shown,  in  Appendix  A, 
that  the  amplitudes  A  and  B  are  related  to  the  amplitudes  F  and  G  by  a  matrix  involving  k,  ic,  and 
the  barrier  width  a.  This  matrix  is  called  the  overall  transfer  matrix,  M.   *:and  k,  of  course, 
depend  directly  on  the  particle  mass  aw,  the  total  particle  energy  £",  and  the  barrier  height  V0,  as 
shown  previously.  In  short,  there  are  four  physical  parameters  which  determine  M:  particle  mass, 
energy,  barrier  height,  and  barrier  width.  The  form  of  the  overall  transfer  matrix  for  the  case  of  a 
simple  rectangular  potential  is  given  in  section  A  of  Appendix  A.  The  overall  transfer  matrix  has 
some  interesting  symmetry  properties;  these  are  developed  in  Appendix  A.  The  determinant  of 
the  overall  transfer  matrix  in  also  discussed  in  Appendix  A. 


ffl.  THE  NUMERICAL  METHODS 

A.  OVERVIEW 

Quantum-mechanical  transmission  through  a  potential  barrier  is  of  interest.  Current 
microprocessors  have  a  device  scale  of  0.3  u.m.  The  room-temperature  de  Broglie  wavelength  of 
a  conduction  electron  in  silicon  is  approximately  0.025  |im.  Future  device  scale  reductions,  by 
merely  a  factor  often,  will  therefore  introduce  significant  quantum  behavior.  If  the  current  pace 
of  microprocessor  scale  reduction  continues,  knowledge  of  the  quantum-mechanical  transmission 
coefficient  will  be  required  in  the  design  of  electronic  devices  within  the  next  15  years. 
(Semiconductor  Industry  Association,  1994)  The  physical  meaning  of  the  transmission  coefficient 
is  described  in  full  by  Eisberg  and  Resnick.  (Eisberg  &  Resnick,  1985,  pp.  193-198)  Essentially, 
the  transmission  coefficient  is  a  measure  of  the  probability  of  barrier  penetration  by  a  quantum 
particle,  such  as  an  electron. 

In  the  preceding  chapter,  we  discussed  the  wave  function  for  a  free  particle.  A  conduction 
electron  in  a  semiconductor,  however,  is  not  a  free  particle:  it  experiences  the  spatially  periodic 
potential  energy  environment  of  the  crystalline  lattice.  Here  we  invoke  the  effective-mass 
approximation,  which,  for  electron  energies  near  the  band  energy  extrema,  permits  us  to  treat  the 
motion  of  the  electron  in  the  solid  as  if  it  were  a  free  particle  having  an  effective  mass,  m*.  The 
value  of  the  effective  mass  is  a  material-specific  parameter.  For  most  devices  of  technological 
interest,  the  electron  energies  remain  near  the  band  extrema  and  the  effective-mass  approximation 
is  valid.  A  complete  description  of  the  effective  mass  is  given  by  Eisberg  and  Resnick.  (Eisberg 
&  Resnick,  1985,  pp.  460-464) 


The  transmission  coefficient  is  defined  as 


T  = 


(10) 


where  the  coefficients  F  and  A  are  the  magnitudes  of  the  transmitted  wave  function  and  the 
incident  wave  function,  respectively.  Let  us  consider  a  barrier  which  extends  from  x=Xl  to  x=Xr. 
The  potential  outside  the  barrier  is  constant.  The  incident  and  reflected  wave  functions  have  the 
form 


V  reflect)  =  Be~ 


*.*  01) 


for  x<xl,  and  the  transmitted  wave  function  is  of  the  form 

W— „*(*)  =  F*«t      (12) 
for  x>Xr,  where  ki  is  the  wave  number  in  the  incident  medium  (before  the  barrier)  and  k2  is  the 
wave  number  in  the  transmitted  medium  (after  the  barrier).  The  wave  number  is  a  complex- 
valued  quantity  defined  by  k=ix,  and  at  is  given  by  equation  9.  Substituting  the  effective  mass, 
m*,  for  the  mass  m,  results  in  the  equation 


2m'(E-V) 

where  V\s  the  constant  value  of  the  potential  energy. 

The  two  numerical  methods  for  finding  the  transmission  coefficient  described  herein  both 
use  solutions  based  upon  piecewise-constant  potentials.  The  first  numerical  method  is  the  transfer 
matrix  method,  which  starts  by  breaking  up  the  actual  potential  into  an  approximate  staircase 
potential,  which  is  piecewise  constant.  A  two-by-two  transfer  matrix  (N  matrix)  is  then  calculated 
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at  each  interface,  and  an  overall  transfer  matrix  (M  matrix)  for  the  entire  potential  is  found  by 
matrix  multiplication  of  the  individual  interfaces'  matrices.  The  transmission  coefficient  for  the 
entire  potential  can  then  be  found.  A  more  detailed  description  of  this  procedure  will  be  given 
later. 

The  second  method  is  a  backward-recurrence  algorithm,  suggested  by  Luscombe 
(Luscombe,  1992,  pp.  1-20),  which  greatly  increases  numerical  stability  by  eliminating  the  error 
associated  with  the  forward  propagation  of  a  solution  involving  both  decaying  and  growing 
exponential  components.  A  complete  description  of  this  technique  also  appears  later. 

This  thesis  emphasizes  the  applications  and  numerical  accuracy  of  these  two  techniques. 
It  will  be  shown  that  the  transfer  matrix  technique  has  a  convenient  built-in  test  for  numerical 
instability.  The  mechanism  of  the  numerical  instability  will  be  discussed  in  detail,  as  will  the 
conditions  which  bring  it  about. 

The  backward-recurrence  method  has  excellent  numerical  stability,  and  also  runs  faster.  It 
does  not  have  a  comparable  built-in  method  to  test  for  accuracy,  however.  The  results  of 
applying  the  method  to  complex  potentials  will  be  discussed. 

There  exists  a  closed-form  analytical  solution  to  the  simple  tunneling  problem,  for  the  case 
of  one  square  barrier  in  a  region  of  otherwise  constant  potential.  (Singh,  1997,  pp.  131-135)  To 
test  the  accuracy  of  both  numerical  methods,  each  technique's  accuracy  in  giving  a  numerical 
solution  to  this  known  problem  was  checked. 

All  computer  programs  used  in  this  thesis  were  written  in  ANSI-standard  C.  This 
language  allows  the  specification  of  either  single-  or  double-precision  floating  point  numbers,  and 
thus  provides  the  ability  to  investigate  numerical  and  algorithmic  stability  at  two  different 
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precision  levels.   The  MATLAB  programming  language  (which  has  built-in  matrix  functionality) 
was  not  used  in  this  work,  primarily  because  it  lacks  this  capability,  but  also  due  to  its  slower 
operation  in  situations  requiring  loops,  which  prove  to  be  unavoidable  in  finding  the  transmission 
coefficient  at  multiple  values  of  incident  particle  energy. 
B.  TRANSFER  MATRIX  METHOD 

The  interface  transfer  matrices  were  the  subject  of  a  particularly  concise  treatment  by 
Singh,  which  I  paraphrase  here:  (Singh,  1997,  pp.  148-149) 

The  transfer  matrix  method  assumes  a  piecewise-constant  potential  approximation  to  the 
potential  of  interest.  Therefore,  the  wave  function  has  the  form 

y/(x)  =  Ae,kx  +  Be',kx       (14) 
in  each  region  of  constant  potential,  as  shown  earlier.  The  incident  wave  is  assumed  to  be 
travelling  in  the  positive  x-direction.  Each  region  will  have  its  own  set  of  coefficients  A  and  B  and 
wave  numbers  k.  Assuming  that  the  electron's  effective  mass  is  constant  in  all  regions,  and  again 
applying  the  boundary  conditions  at  the  interface  between  regions,  here  situated  at  x=x0. 


V(0  =  V(0 

dx  x°       dx  *° 


(15) 


we  see  that  the  coefficients  A  and  B  can  be  related  to  each  other  by  the  matrix  equation 


B, 


=  N 


A2 
B, 


(16) 


where  A\  and  B\  are  the  coefficients  of  the  wave  function  for  x  x0,  and  A?  and  B2  are  the 
corresponding  coefficients  for  x  x0.  The  matrix  N  is  called  the  interface  transfer  matrix.  This 
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two-by-two  matrix  can  be  shown,  in  the  general  case  of  a  potential  with  interfaces  at  points  x, 
separating  regions  of  potential  Vi+i,  to  have  components 


N,(U) 


K(12)  = 


1+1 


1  + 

v      K 


,teX*w-*i! 


1     k,+x  \c(-'x')[k"l+K) 


\ 


V 


1     ^'+1  U"' )<*»+■ +*'> 


(17) 


•■!  +  l 


1  + 


,(-*;-)(  *w-*iJ 


The  matrix  product  of  all  of  the  N  matrices  yields  the  total  transfer  matrix  for  the  potential,  M. 
We  use  the  total  transfer  matrix  to  calculate  the  overall  transmission  coefficient,  by  the  relation 

1 


T  = 


|M(1,1)| 


(18) 


This  method  is  robust  for  simple  potentials,  but  runs  into  difficulties  if  the  barrier  height,  barrier 

width,  or  number  of  interfaces  is  too  great.  This  will  be  discussed  in  greater  detail  later.  The 

program  shows  perfect  agreement  with  the  analytic  form  of  the  transmission  coefficient  for  one 

square  barrier. 

C.  BACKWARD-RECURRENCE  METHOD 

The  backward-recurrence  method  used  in  this  thesis  is  that  described  by  Luscombe. 

(Luscombe,  1992,  pp.  1-20)  As  he  states, 

relevant  device  variables  must  ...  be  specified  within  close  tolerances,  and  clearly 
the  most  efficient  means  of  developing  a  realistic  [microprocessor]  design  is 
through  computer  modeling  based  on  fundamental  physical  laws,  e.g.  the  Poisson 
and  Schrodinger  equations    By  providing  the  tools  to  explore,  in  detail,  and 
before  fabrication,  the  effects  of  varying  ...  [device]  parameters  ...,  modeling 
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optimizes  the  development  cycle  and  provides  a  cost-effective  method  for 
validating  and  refining  device  concepts.  (Luscombe,  1992,  p. 2) 


Before  discussing  the  backward-recurrence  method  in  detail,  we  note  the  following 
observations.  The  Schrodinger  equation  is  a  second-order  differential  equation,  and,  as  such,  has 
two  linearly  independent  solutions.  In  the  classically  forbidden  potential  regions,  the  physical 
solution  for  the  wave  function  ifr(x)  is  strictly  decreasing  as  a  function  of  distance.  (Since  E<  Fin 
the  classically  forbidden  regions,  the  second  derivative  of  ^is  strictly  positive.)  Now,  as  is  simple 
to  check,  for  example  by  examining  the  Wronskian  relation  associated  with  the  Schrodinger 
equation,  if  one  solution  is  monotonically  decreasing,  the  other,  linearly  independent  solution  is 
monotonically  increasing.  This  second,  linearly  independent  solution  is  therefore  not  physical.  In 
numerical  methods  that  propagate,  through  the  use  of  iteration,  both  the  growing  and  decaying 
solutions  on  equal  footing  (as  in  the  transfer  matrix  method),  the  slightest  round-off  error  will 
trigger  the  growth  of  the  unwanted,  growing  solution  in  the  classically  forbidden  potential 
regions.  By  employing  backward  iteration,  however,  the  physical  solution  becomes  dominant, 
increasing  resolution  in  the  direction  of  backward  iteration. 

The  transfer  matrix  method,  while  accurate  in  limited  applications,  is  numerically  unstable, 
as  stated  above.  This  instability  is  implicit  in  the  N  matrices  (interface  matrices)  described  earlier, 
(discussions  with  James  Luscombe,  May- June  1997) 

The  backward-recurrence  method  works  as  follows.  We  shall  assume  the  electron's 
effective  mass  to  be  constant  throughout,  and  equal  to 

m'=0.067me;       (19) 
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we  therefore  need  not  change  it  from  region  to  region  (me,  here,  is  the  electron  mass).  The  only 

thing  that  changes  spatially  is  the  potential,  V.  As  with  any  numerical  solution  to  a  differential 

equation,  we  employ  the  finite  difference  approximation.  We  sample  the  wave  function  at  discrete 

points  x„=nA,  where  A  is  the  step  size,  and  is  chosen  to  be  sufficiently  small  to  accurately  resolve 

the  potential.  We  begin  with  the  incremental  Taylor  series  expansions 

A2 
y/(x  +  A)  =  y/(x)  +  A  y/'  ( x)  H y/"(x)  +  (higher  -  order  terms) 

2  (20) 

A2 
y/(x  -  A)  =  y/(x)  -  Ay/'(x)  -\ yf"{x)  +  (higher  -  order  terms). 


An  expression  for  y/", 


ys"(x)  =  ^ )—^-2 1 ^,       (21) 


can  then  be  derived.  The  algorithm  progresses  from  right  to  left,  through  the  perturbation 
described  below,  in  spatial  steps  of  size  A.  The  solution  at  the  11th  step  may  be  written  as 

y/=^,+yr.-2yv    (22) 

The  time-independent  Schrodinger  equation, 

-^Vn+Vy/  =  Eyf,      (23) 
2/w 

is  equivalent  to  a  three-term  recurrence  relation, 

V^.+Wi+^Vn^       (24) 


where  bn  is  given  by 


^2m*An 

h2 


(E-V„)-2.       (25) 
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Finally,  using  the  definition 

r„  =  -^,      (26) 

the  above  three-term  recurrence  relation  can  be  rewritten  as 

r»=rr— '    (27) 

which  is  easily  recognizable  as  a  two-term  backward-recurrence  relation.  It  is  this  relation  which 
is  the  heart  of  the  backward-recurrence  method. 

The  potential  being  analyzed  must  consist  of  an  irregularity  of  finite  spatial  dimensions, 
beyond  which  the  potential  is  constant.  The  irregularity  in  the  potential  may  be  of  arbitrary  height 
(eV),  curvature,  and  width  (nm).  There  must  also  exist  two  regions  of  constant  potential,  which 
we  designate  as  right-hand  and  left-hand  regions.  Each  has  an  associated  wave  number,  kR  and  kL. 
The  constant  potential  energies  in  these  two  regions  need  not  be  the  same,  necessarily,  but  for 
simplicity  they  are  assumed  to  be  identical  in  this  thesis. 

Defining  r0  to  be  the  value  of  r  at  the  edge  of  the  left-hand  region,  it  is  easily  seen  that 

„=_^      ^[(0)(A)]  _         1+j 

•     r-.     (4(-l)(A)]     e-**+Re*>"      (28) 

where  t/r„= i/r(x)  =e'kLx  +  e'     and  x=nA,  and  that  the  reflection  coefficient  R  is  given  by 

re~*iA  -1 
R  =  — rr-         (29) 

o 

Similar  logic  can  be  applied  to  the  right-hand  region  of  constant  potential,  by  writing  rN+i  as 
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¥nu  _  Te 


:kRA(N+\) 


A/  +  1 


¥, 


Te 


ikRA(N) 


=  e'k«&  =  r, 


N  ' 


(30) 


ik„x 


where  if/ n  =  y/(x)  =  Te  rX  . 


We  know  rN+!  if  we  know  kR  and  the  step  size,  A.  The  expression  for  kR  is 


**  = 


2m 


-U-VR),       (31) 


where  VR  is  the  constant  potential  in  the  right-hand  region.  An  analogous  expression  exists  for  kL, 


h  = 


(irn^ 


^    J 


(E-VL).     (32) 


We  have  derived  one  expression  for  R,  Equation  29.   T\s  related  to  R  by  the  expression 

|7f  =  !J5M(l-)*fV       (33) 

cm  t    A  \  / 


sin/:RA 


where  T  and  R  are  the  overall  transmission  and  reflection  coefficients,  respectively.  (Luscombe, 
1992,  p.  7)  Equation  33  may  be  rewritten  as 


\jf  = 


rsinkLA}(    (4  sin  kLA)  Im(r0) 
smkRA 


l-2Re(^A)  +  |r0| 


(34) 


which  the  program  uses  to  calculate  T  once  it  has  found  the  value  of  ra. 
D.  ANALYTIC  SOLUTION 

The  analytic  solution  to  the  problem  of  one  square  barrier,  used  to  test  the  transfer  matrix 
and  backward-recurrence  methods,  is  simple  to  derive.  (Eisberg  &  Resnick,  1985,  pp.  199-201) 
(Singh,  1997,  pp.  131-134)  Consider  a  square  barrier  of  height  V0  and  width  a,  extending  from 
x=0  to  x=a.  The  potential  everywhere  outside  the  barrier  is  zero.  Define  the  region  to  the  left 
and  right  of  the  barrier  to  have  wave  number  &/,  and  the  barrier  to  have  associated  wave  number 
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kn.  These  wave  numbers  are  given  by 


2m 


rlE-o) 


'2m*^ 


(35) 


(E-K). 


V  "    ) 


As  stated  earlier,  the  solution  to  the  Schrodinger  wave  equation  in  a  region  of  constant  potential 

has  the  form  y/{x)  =  Aelkx  +  Be~,kx  .  Call  the  region  in  which  x <0  region  I,  the  region  in  which 

0<x<a  region  II,  and  the  region  in  which  x   a  region  III.  The  wave  functions  in  these  regions  are 

y,  =Ae'k'x+Be~,k'x 

y/u  =Ce*»x+De-*°x(36) 

¥m=Feik>x, 

assuming  that  there  is  only  a  right-moving  wave  in  region  III.  Applying  the  boundary  conditions 
on  the  wave  functions  given  by  Equations  15,  at  x=0  and  x=a,  and  performing  algebraic 
simplification,  one  arrives  at  the  expression 


T2  = 


T"/C  t  f\  f  7  t- 


ia(ku-k, 


(*,+*„)   -fa-tn)*"* 


(37) 


which  is  the  expression  plotted  by  the  programs  to  check  the  accuracy  of  the  two  numerical 
methods. 
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IV.  THE  PROGRAMS 

A.  TRANSFER  MATRIX  PROGRAM 

A  copy  of  the  transfer  matrix  program,  method  1  c,  is  attached  as  Appendix  C,  Section  A. 
The  program  is  written  in  ANSI  C. 

Since  ANSI  C  contains  no  built-in  complex  numbers  capability,  a  global  variable  structure 
called  complex  is  defined,  consisting  of  two  double-precision  floating  point  numbers  representing 
the  real  and  imaginary  parts  of  the  complex  number.  Also,  since  the  potential,  V,  and  wave 
number,  k,  as  functions  of  x,  are  required  in  multiple  functions,  they  are  defined  globally.  All 
other  variables  have  local  scope. 

The  functions  addc,  subc,  mule,  and  dive  have  been  written  to  perform  the  four  basic 
mathematical  functions  on  complex  arguments.  They  all  take  two  arguments  of  type  complex, 
and  return  type  complex.  The  function  expc  computes  ex,  where  x  is  a  complex  argument,  and 
returns  type  complex.  The  function  absc,  given  a  complex  argument,  returns  the  (real)  magnitude 
of  the  complex  number,  as  a  double-precision  floating-point  number.  The  function 
determinantcomplex  takes  four  arguments  of  type  complex,  which  represent  the  elements  of  a 
two-by-two  overall  transfer  matrix,  and  returns  the  (complex)  determinant  of  the  matrix. 

Nl  1  complex,  N12complex,  N21  complex,  and  N22complex  calculate  the  values  of  the 
elements  of  the  interface  transfer  matrix,  using  Equation  17.  They  take  four  arguments  of  type 
complex,  representing  the  wave  vector  left  and  right  of  the  interface,  their  ratio,  and  the  value  of  x 
at  the  interface. 
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Finally,  the  function  makeV  initializes  the  potential  energy  array,  V.  In  this  program,  the 
potential  energy  profile  is  a  series  of  square  barriers  of  height  Vo,  width  BARWIDTH,  and 
separated  by  BARWIDTH.  BARRIERS  represents  the  number  of  barriers  present.  ERIGHT, 
ELEFT,  and  EPTS  represent  the  highest  and  lowest  values  of  incident  particle  energy,  and  the 
number  of  energy  values  used,  respectively    Vo,  BARWIDTH,  BARRIERS,  ERIGHT,  ELEFT, 
and  EPTS  are  all  predefined  constants,  listed  as  #define  statements. 

Three  other  physical  constants  are  used.  One  is  the  effective  mass,  defined  as  a  ratio, 

m 

— ,  equal  to  0.067,  and  represented  by  EFFMASS.  This  is  an  appropriate  value  of  the  effective 

me 

h2  2 

mass  for  GaAs  semiconductors.  The  value  of ,  in  units  of  eVnm  ,  is  defined  by  H20VER2M 

2m 

e 

as  0.0381,  and  is  constant  for  all  materials.  Finally,  the  program  calculates  a  constant,  C,  which  is 

,       2m      „    L  .      ,  ,  J         EFFMASS 

equal  to  — =-  .  C,  obviously,  may  be  expressed  as 


h-  r  H20VER2M 

Function  main  generates  the  program's  output,  and  controls  the  execution  of  the  transfer 
matrix  method  algorithm.   The  printf("%.  12f\t%.  12f\n",E,T),  statement  lists  the  values  of  incident 
particle  energy  and  transmission  coefficient  in  two  columns,  to  the  screen.  To  generate  the  graphs 
found  herein,  the  program's  output  has  been  saved  in  a  file  using  indirection,  and  plotted  using  the 
commercially  available  program  Spyglass  Plot. 

The  basic  execution  of  main  is  as  follows.  A  value  of  incident  particle  energy  is  assigned 
by  the  for(ec=0.0;ec<epts;ec+=l  .0)  statement.  This  loop  goes  through  all  desired  values  of 
energy,  from  ERIGHT  down  to  ELEFT,  dividing  the  energies  into  EPTS  steps.  As  EPTS  is  an 
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integer,  the  variable  epts  is  introduced  to  promote  it  to  a  double-precision  floating-point  number 

in  order  to  be  type  compatible  with  the  variable  ec.  This  also  prevents  errors,  which  might  be 

caused  by  inadvertent  integer  division. 

Once  a  value  of  £  has  been  defined  by  the  E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); 

statement,  the  variables  Ml  /,  Ml  2,  M21,  and  M22,  all  of  type  complex,  are  initialized.  These 

variables  are  the  matrix  elements  of  the  overall  transfer  matrix.   Since  the  overall  transfer  matrix 

will  be  calculated  by  matrix  multiplication  of  all  intervening  interface  transfer  matrices,  it  is 

initialized  as  a  two-by-two  identity  matrix  by  the 

oldMl  1  .imag=oldM  1 2.imag=oldM2 1  .imag=oldM22.imag=0.0; 
oldMl  l.real=oldM22.real=1.0; 
oldM12.real=oldM21.real=0.0; 

statements. 

The  next  step  is  to  initialize  the  wave  vector  array,  k,  at  the  current  value  of  E.  The 
for(xc=0;xc<XPTS;xc++)  loop  handles  this.  At  each  point,  the  wave  vector  will  be  either  purely 
real  or  purely  imaginary,  and  is  given  by  Equation  1 3 . 

The  for(xc=l;xc<XPTS;xc++)  loop,  which  follows,  goes  through  the  potential  profile  and 
identifies  interfaces;  that  is,  it  finds  values  of  x  at  which  k  changes.  When  an  interface  is  found, 
the  if((k[xc].real!=k[xc-l] .real)||(k[xc].imag!=k[xc-l].imag))  statement  executes,  causing  the 
interface  transfer  matrix  to  be  calculated  and  matrix  multiplied  by  the  "running  total"  overall 
transfer  matrix,  given  by  elements  oldMll,  oldM12,  oldM21,  and  oldM22,  all  of  type  complex. 
Upon  completion  of  this  if  statement,  the  program  checks  for  numerical  instability. 

The  Mdet=determinantcomplex(newMl  l,newM12,newM21,newM22);  statement  checks 
the  determinant  of  the  calculated  overall  transfer  matrix.  The  subsequent 
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if((Mdet.real>l.  000  l)||(Mdet.real<0. 9999))  statement  checks  to  see  if  the  determinant  is  one.  If  it 
is  not,  within  four-digit  precision,  then  the  program  halts  execution  and  prints  the  listed  error 
message. 

The  transmission  coefficient,  T,  is  then  found  using  Equation  1 8,  by  the 
T=sqrt(1.0/(absc(newMl  l)*absc(newMl  1)));  statement.  It  is  then  printed  out  in  two  columns,  as 
mentioned  before.   Spyglass  Plot  then  produces  graphs  from  the  output. 
B.  BACKWARD-RECURRENCE  PROGRAM 

A  copy  of  method2.c,  which  employs  the  backward-recurrence  method,  also  is  included 
in  Appendix  C,  Section  B.  The  functions  addc,  subc,  mule,  dive,  expc,  and  absc  are  re-used  here, 
to  perform  the  basic  operations  on  complex  numbers.  The  predefined  constants  are  also  the  same 
as  those  used  in  method  1  c,  with  the  exception  that  VL,  fM,  and  VR  have  been  added  to  allow  the 
specification  of  different  potential  in  the  left-hand  and  right-hand  constant-potential  regions,  and 
between  the  potential  barriers.  The  program  operates  as  follows: 

Variables  r,  ro,  rnplusJ,  and  bnfXPTSJ,  of  type  complex,  hold  the  values  required  in  the 
backward-recurrence  relation  given  by  Equation  27.  Function  makeV,  as  before,  creates  a 
potential  profile  consisting  of  square  barriers.  The  barriers  have  height  Vo,  and  separation  and 
width  given  by  BARW1DTH .  Other  versions  of  this  program  have  been  written  which  calculate 
the  transmission  coefficient  through  different  potentials.  Only  makeV  need  be  changed  to 
accomplish  this. 

The  for(ec=0.0;ec<epts;ec+=1.0)  loop  counts  through  the  values  of  incident  particle 
energy,  from  ERIGHT  down  to  ELEFT,  using  the  E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); 
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statement.  At  each  value  of  incident  particle  energy,  the  for(xc=0;xc<XPTS;xc++)  loop  initializes 
the  bn  array,  bn  is  given  by  Equation  25,  which  is 

bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 

bn[xc].imag=0.0, 

in  the  C  language. 

kl  and  kr  represent  the  wave  numbers  in  the  left-hand  and  right-hand  regions  of  constant 

potential,  respectively,  and  are  given  by  Equations  3 1  and  32  and  by 

kl=sqrt(C*(E-VL)); 
kr=sqrt(C*(E-VR)); 

Equation  30  gives  rN+!,  the  value  of  r  in  the  right-hand  region  of  constant  potential.  The  variable 

rnplusl  represents  rN+,  in  the  program,  and  it  is  initialized  by  the 

rnplusl  real=cos(kr*delta); 
rnplusl  imag=sin(kr* delta); 

statements. 

The  for(xc=(XPTS-BARPTS);xc>=BARPTS;xc--)  loop,  while  deceptively  short,  actually 

performs  all  of  the  backward-recurrence  steps,  using  Equation  27,  and  updating  the  value  of 

rnplusl  before  each  iteration.  Finally,  having  computed  ro,  the  program  calculates  the  value  of 

the  transmission  coefficient,  T,  using  Equation  34.  This  is  done  by  the 

temp.real=0.0; 
temp .  imag=kl  *  delta; 
temp2=mulc(ro,expc(temp)); 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 
2.0*temp2.real+absc(ro)*absc(ro))); 
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statements.  The  values  of  E  and  Tare  then  printed  out,  in  columns,  as  before,  and  the  results 
plotted  using  Spyglass  Plot. 
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V.  PERFORMANCE  AND  NUMERICAL  STABILITY 

A.  TRANSFER  MATRLX  METHOD 

The  transfer  matrix  method,  as  mentioned,  contains  a  testable  quantity  which  allows  one 
to  determine  when  its  output  has  begun  to  be  unreliable.  It  is  shown  in  Section  B  of  Appendix  A 
that  the  determinant  of  the  overall  transfer  matrix,  M,  for  a  system  must  equal  one.  This  has  been 
used  to  find  the  point  at  which  the  program's  output  begins  to  deteriorate. 

By  testing  the  determinant  in  this  way,  one  also  may  find  the  specific  value  of  incident 
particle  energy  at  which  the  transfer  matrix  method  breaks  down.  Extensive  results  of  these  tests 
are  attached  as  Appendix  B,  for  a  simple  test  potential  composed  of  a  series  of  square  barriers, 
whose  number,  height,  and  width/separation  may  be  changed  until  the  program  detects  numerical 
instability.  From  these  test  cases,  it  is  clear  that  the  transfer  matrix  method  can  be  numerically 
unstable  in  many  applications.  Examples  of  this  instability  will  be  given  shortly. 

In  Appendix  B,  Section  A  is  the  output  of  the  transfer  matrix  method  program, 
methodic,  when  single-precision  floating-point  numbers  are  used.  These  numbers  have  the 
equivalent  of  six  decimal  places  of  precision.  Consider  the  case  in  Appendix  B,  Section  A,  Part  1 
(B.A.  1 .) ,  of  one  five  nanometer  (nm)  barrier  which  is  0.2  electron  volts  (eV)  tall.  In  this  case,  for 
single-precision  numbers,  the  method  breaks  down  at  an  incident  particle  energy  of  0.010702  eV. 
method l.c  analyzes  particle  energies  beginning  with  the  largest  value,  so  for  any  particle  energy 
less  than  this  value,  the  transfer  matrix  method  returns  inaccurate  results.  Compare  this  result  to 
that  in  B.B.  1.,  which  is  the  same  physical  case,  analyzed  using  double-precision  floating-point 
numbers,  which  carry  12  digits  of  precision.  No  numerical  instability  occurs  with  the  double- 
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precision  numbers;  this  is  not  surprising,  since  more  precision  is  available  in  the  calculations  used 
in  the  program.  There  is  no  number  of  digits  of  precision,  however,  which  is  sufficient  to  entirely 
prevent  numerical  instability  of  the  transfer  matrix  method.  Note  in  B.B.  1.,  for  instance,  that 
inaccuracy  occurs  when  the  height  of  a  single,  5  nm  barrier  is  increased  to  5  eV. 

In  B.A.2.  and  B.B. 2.,  the  barrier  width  has  been  varied.  Note  that  for  both  single-  and 
double-precision  floating-point  numbers,  there  exists  a  value  of  barrier  width  which  causes 
numerical  instability  of  the  transfer  matrix  method.   Similarly,  in  B.A.3.  and  B.B. 3.,  this  is 
observed  when  the  number  of  square  potential  energy  barriers  is  increased.  In  each  case,  the 
value  of  incident  particle  energy  at  which  instability  occurs  is  listed.  Observe  that,  as  the 
parameters  of  the  calculation  are  pushed  beyond  the  point  at  which  instability  just  starts,  the 
minimum  incident  particle  energy  rises.  This  effect  can  be  seen  for  either  single-  or  double- 
precision  floating-point  numbers,  throughout  Appendix  B. 

The  question  arises:  how  severe  is  this  effect  in  the  transfer  matrix  method?  At  this  point, 
methodic  is  applied  to  a  situation  with  a  known,  analytic  solution,  in  order  to  find  out  the  answer 
to  this  question.  The  following  paragraphs  will  refer  to  various  figures,  all  of  which  have  been 
collated  in  Appendix  D. 

The  test  case  chosen  is  that  of  one  square  potential  energy  barrier,  of  width  5.0  nm  and 
height  0.23  eV.  In  B.A.2.,  it  has  been  shown  that  the  transfer  matrix  method  breaks  down  at 
incident  particle  energy  0.013768  eV,  for  single-precision  numbers.  It  does  not  fail  if  double- 
precision  numbers  are  used.  Figure  1  shows  the  output  of  the  program  using  double-precision 
numbers.  This  figure  is  a  familiar  sight  in  introductory  quantum  physics  texts.  (Eisberg  & 
Resnick,  1985,  p.  202)  (Singh,  1997,  p.  135)  There  is  a  positive  transmission  coefficient  for 
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incident  particle  energies  less  than  the  barrier  height,  indicating  that  quantum  tunneling  will  occur. 
There  is  also  a  region  of  non-unity  transmission  coefficient  at  values  above  the  barrier  height, 
indicating  the  presence  of  quantum  scattering. 

In  Figure  2,  the  transfer  matrix  method  solution  minus  the  value  of  the  transmission 
coefficient  calculated  by  the  analytic  solution  for  one  square  barrier  has  been  plotted.  Note  that 
the  transfer  matrix  solution  is  exact.  The  same  figure  resulted  from  running  the  program  with 
single-precision  floating-point  numbers.  It  is  apparent  that,  when  the  determinant  of  the  overall 
transfer  matrix  is  1.0001  instead  of  1.0000,  significant  degradation  of  the  numerical  solution  is  not 
observed. 

This  is  an  interesting  dilemma,  since  analytic  solutions  for  the  transmission  coefficient  are 
so  rare  in  practical  quantum  tunneling  problems.  When  is  the  output  of  the  transfer  matrix 
method  useful?  What  are  the  limitations  of  the  method?  Only  qualitative  answers  are  available  to 
the  first  question;  the  method  is  useful  when  it  is  numerically  stable.  Instability  may  be  detected 
by  means  of  the  determinant  of  the  overall  transfer  matrix.  The  limitations  of  the  transfer  matrix 
method,  however,  are  tested  in  this  thesis  by  five  test  potentials:  a  single  square  barrier,  sequences 
of  two,  three,  and  five  parabolic  barriers,  and  a  resonant-tunneling  diode  (RTD)  potential. 

The  results,  in  the  case  of  the  RTD  potential,  clearly  demonstrate  the  inadequacy  of  the 
transfer  matrix  method  in  dealing  with  elaborate  potential  profiles.  Appendix  C,  Section  E 
describes  ml  para,  c,  which  tests  the  transfer  matrix  method  against  potential  energy  profiles 
composed  of  series  of  parabolic  barriers,  rather  than  square  barriers.  An  example  of  a  potential 
used  in  this  program  is  given  in  Figure  3.  In  this  figure,  one  can  clearly  see  the  difficulty  posed  by 
potentials  consisting  of  smooth  curves:  the  potential  must  be  carved  into  sub-bins,  each  of  which 
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is  piecewise-constant.  The  partitioning  must  be  such  that  the  character  of  the  potential  in 
preserved,  which,  for  steep  slopes  in  the  potential,  may  require  small  step  sizes.  Recall  that  the 
number  of  bins  (barriers)  tends  to  drive  the  method  into  instability,  as  shown  previously. 

Figures  4,  5,  and  6  are  the  output  of  ml  para,  c  for  the  cases  of  2,  3,  and  5  parabolic 
barriers.  In  each  of  these  cases,  the  "safety  check"  of  the  determinant  of  the  overall  transfer 
matrix  has  been  performed  for  each  value  of  energy,  and  program  execution  halted  when  the 
determinant  differs  from  1 .0  by  more  than  0.0001 .  In  the  case  shown  in  Figure  4,  instability  did 
not  occur.  In  Figure  5,  it  occurred  after  952  out  of  1000  energy  points  had  been  analyzed.  In 
Figure  6,  it  happened  after  only  908.  Each  graph,  however,  shows  only  the  values  of  T  which 
were  calculated  before  instability  occurred.  These  values  are  therefore  known  to  be  accurate. 

It  has  been  established  that  the  output  of  the  transfer  matrix  method  is  correct,  even  when 
the  determinant  of  the  overall  transfer  matrix  differs  slightly  from  one.  When  the  discrepancy  in 
the  determinant  is  less  than  0.0001,  therefore,  the  program's  output  is  accurate.  These  figures 
will  be  used  as  a  basis  for  evaluating  the  performance  of  the  backward-recurrence  method,  applied 
to  the  same  potentials,  in  Section  B  of  this  chapter. 

The  results  seen  in  Figure  4  seem  reasonable.  There  is  a  peak  transmission  coefficient  of 
1.0  when  the  energy  equals  the  height  of  the  barriers,  which  agrees  with  observations  based  on 
systems  of  square  barriers.   Scattering  occurs  for  incident  energies  above  the  barrier  energy,  as 
expected.   Since  there  are  two  barriers,  one  expects  a  resonant-tunneling  energy  to  exist,  and 
there  is,  in  fact,  a  peak  in  the  T  versus  E  curve  for  incident  energies  below  the  barrier  energy.  The 
performance  of  the  transfer  matrix  method  in  this  case  is  acceptable. 
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Figure  5,  for  the  case  of  three  parabolic  barriers,  demonstrates  that  the  technique  is  fairly 
accurate,  as  well.  There  is  a  scattering  region  above  the  barrier  energy,  as  expected,  and  there  is 
also  a  resonant-tunneling  energy.  In  addition,  we  see  the  expected  twinning  of  transmission 
resonances  due  to  the  presence  of  two  energy  wells  between  the  three  barriers.  However,  as  noted 
above,  the  method  breaks  down  at  low  energies. 

In  Figure  6,  similar  performance  occurs.  This  figure  is  plotted  with  a  logarithmic  vertical 
axis,  to  better  show  the  range  of  values.  Again,  twinning  of  transmission  resonances  appears, 
with  four  resonances  this  time.  This  is  due  to  the  four  energy  wells  that  exist  between  the  five 
energy  barriers.  Performance  degrades  at  low  energies  once  more,  and  instability  sets  in  at  lower 
energies  this  time,  as  anticipated. 

Finally,  the  transfer  matrix  method  is  applied  to  a  truly  difficult  (and  physically  relevant) 
potential:  the  resonant  tunneling  diode  (RTD)  potential.  This  potential  is  shown  in  Figure  7,  for 
an  AlAs/InGaAs/InAs  RTD.  (Luscombe,  1992,  p.  4)  Note  the  steeply-sloping  sides  of  the  central 
well,  in  combination  with  the  smooth  curve  on  which  the  central  well  rests.  These  features  make 
it  extremely  difficult  to  approximate  this  potential  by  a  piecewise-constant  model. 

Figure  8  shows  the  results  obtained  from  the  transfer  matrix  method  in  this  case.  The 
program's  output,  even  using  double-precision  numbers,  is  so  inaccurate  that  the  transmission 
coefficient  is  not  calculated  between  0.75  and  0.87  eV,  and  the  method  breaks  down  completely 
below  about  0.55  eV.  Only  the  energies  of  the  two  peaks  are  correctly  predicted,  and  the  method 
completely  misses  a  third,  low-energy  peak,  as  described  in  Section  B  of  this  chapter. 

In  defense  of  the  method,  it  must  be  said  that  for  those  potentials  which  are  simple  enough 
for  it  to  handle,  it  performs  perfectly.  However,  the  transfer  matrix  method  does  have  the 
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observable  disadvantage  of  considerably  longer  execution  time,  for  the  same  potential,  when 
compared  to  the  backward-recurrence  method.  The  observed  time  difference  has  been  as  much  as 
several  minutes,  on  a  Sun  SparcStation  5  workstation,  depending  on  the  specific  problem. 

The  method  also  proves  to  be  able  to  operate  accurately  even  when  the  determinant  of  the 
overall  transfer  matrix  differs  slightly  from  one,  but  it  rapidly  becomes  unstable  for  determinant 
values  which  differ  from  one  by  more  than  0.001 .  Such  differences  may  easily  occur  if  the 
number  of  interfaces  is  greater  than  ten.  Difficulties  arise  when  the  program  is  required  to 
continue  to  handle  many  more,  lower,  incident  particle  energies,  beyond  that  energy  which  first 
caused  the  determinant  not  to  equal  one.  The  incident  energies  are  analyzed  in  decreasing  order, 
so  that  the  method  first  encounters  the  cases  least  likely  to  cause  numerical  instability.  Cases  do 
exist,  however,  for  which  the  transfer  matrix  method  is  grossly  incapable  of  calculating  the 
transmission  coefficient,  and  it  fails  completely.  As  these  cases  are  of  significant  interest  in  the 
development  of  nanoelectronic  devices  (like  the  RTD),  the  transfer  matrix  method  proves  not  to 
be  powerful  enough  for  these  applications. 

B.  BACKWARD-RECURRENCE  METHOD 

The  backward-recurrence  method  has  proven  to  be  considerably  faster,  easier  to  program, 
and  less  prone  to  numerical  instability  than  the  transfer  matrix  method.  In  fact,  it  has  not  been 
demonstrated  to  fail  numerically  at  all,  though  the  fine  detail  of  the  T  versus  E  curves  of  some 
potentials  can  also  be  a  factor,  as  will  be  shown  below,  on  page  32. 
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Since  it  lacks  a  built-in  test  for  numerical  accuracy,  the  best  tests  for  the  backward- 
recurrence  method  are  comparisons  with  known  transmission  coefficient  versus  energy  curves. 
For  example,  again  consider  the  case  of  one  square  barrier. 

In  Figure  9,  the  output  of  the  backward-recurrence  program  has  been  plotted  on  the  same 
axis  as  that  of  the  analytic  expression  for  T  for  one  square  barrier.  Note  the  close  agreement 
between  the  numerical  and  analytic  solutions.  In  Figure  10,  the  difference  between  the  analytic 
result  and  the  backward-recurrence  result  is  shown.  The  agreement  is  perfect,  to  the  limits  of 
precision  of  the  calculation. 

The  only  other  reference  which  one  can  compare  to  the  results  of  the  backward-recurrence 
method  are  known  correct  results  from  the  transfer  matrix  method.  (This  makes  clear  why  the 
transfer  matrix  method  has  been  considered:  it  can  tell  the  user  when  it  has  generated  correct 
results.)  The  T  versus  E  profiles  of  Figures  4  through  6,  generated  using  the  transfer  matrix 
technique,  can  be  used  for  comparison. 

The  case  depicted  in  Figure  1 1  is  the  same  as  that  of  Figure  4,  except  calculated  with  the 
backward-recurrence  technique.  Comparison  of  the  two  graphs  suggests  that  the  agreement  of 
these  two  figures  is  perfect.  From  this  analysis,  it  is  clear  that  the  backward-recurrence  results  and 
the  transfer  matrix  results  agree  to  within  the  calculations'  precision. 

Similar  comparisons  have  been  done  between  the  data  of  Figure  12  and  Figure  5,  and 
between  that  of  Figure  13  and  Figure  6,  for  the  cases  of  three  and  five  parabolic  barriers.  For 
both  of  these  potentials,  the  backward-recurrence  method  exactly  reproduces  the  known-correct 
results  obtained  from  the  transfer  matrix  method.  The  backward-recurrence  method  executes 
these  calculations  in  approximately  one-tenth  the  time  required  by  the  other  algorithm,  as  well. 
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Plot  resolution  impacts  the  quality  of  the  program  output,  even  when  the  individual  data 
pairs  (E,T)  are  correct.  The  rate  of  change  of  T  with  changes  in  E  is  so  sharp  in  the  neighborhood 
of  the  maxima,  that  with  1000  energy  points,  the  peak  itself  is  stepped  over.  Figure  13  contains 
this  error,  which  is  unrelated  to  the  accuracy  of  the  individual  values  of  T.  We  see  in  Figure  13 
that  the  magnitudes  of  T  at  the  four  maxima  located  at  about  0.06  eV  are  not  all  1.0.  They  should 
be,  as  Figure  14,  which  uses  10,000  energy  points  instead  of  1000,  demonstrates. 

As  an  example  of  this  phenomenon,  compare  the  leftmost  maximum  near  0.06  eV  in 
Figure  13,  with  the  same  maximum  in  Figure  14.  In  Figure  13,  the  frequency  at  which  energy 
values  are  sampled  is  insufficient,  and  this  T  maximum  appears  to  be  less  than  1.0.  In  Figure  14, 
the  energy  values  are  sampled  frequently  enough  to  show  the  true  maximum  of  1 .0.  Required 
graphical  resolution  for  plots  of  this  type  may  easily  exceed  the  specified  energy  sample  rate.  This 
highlights  the  importance  of  correct  sample  settings  to  the  proper  use  of  these  numerical 
techniques.  Fortunately,  for  nanoelectronic  devices,  the  number  of  layers,  and  thus  the  number  of 
material  interfaces,  will  be  finite;  this  will  alleviate  some  of  this  problem,  since  it  will  limit  the 
slope  of,  and  the  number  of  maxima  in,  the  T  versus  E  profile  of  the  device. 

Using  backward  recurrence,  the  transmission  coefficient  for  the  RTD  potential  profile  was 
easily  calculated.  These  results  have  been  included  as  Figure  15.  This  figure  was  plotted  with 
10,000  energy  points  after  noting  that  the  maximum  at  about  0.6  eV  had  a  value  of  T  which  was 
less  than  1.0,  when  plotted  using  1000  energy  points.  Note  the  full  coverage,  lack  of  numerical 
instability  at  low  energies,  and  three  energy  peaks.  This  result  is  in  agreement  with  Luscombe's 
original  paper,  taking  into  consideration  differences  in  the  effective  masses  of  the  media. 
(Luscombe,  1992,  p.  4)  The  backward-recurrence  method  proves  itself  to  be  quite  capable  of 
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attacking  the  RTD  potential,  and  it  no  doubt  is  capable  of  analyzing  the  potential  profiles  of  other 
nanoelectronic  devices.   This  method  has  great  promise,  and  wide  future  applicability. 
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VI.  CONCLUSION 

The  backward-recurrence  method  has  much  greater  stability  and  much  faster  speed  of 
execution,  as  well  as  a  much  wider  scope  of  applicability  than  does  the  transfer  matrix  approach. 
Both  methods  reproduce,  with  no  observable  error,  the  transmission  coefficient  (7)  versus 
incident  particle  energy  (E)  curve  for  the  classic  square  potential  barrier,  as  seen  in  Figures  1  and 
2. 

That  the  determinant  of  the  transfer  matrix  equals  one,  means  that  the  results  produced  by 
the  transfer  matrix  method  are  credible.  This  fact  has  been  used  to  calculate  T  versus  E  curves  for 
other  potentials,  such  as  the  parabolic  barrier  potentials  seen  in  Figures  4  through  6.  These  T 
versus  E  curves  found  with  the  transfer  matrix  method  have  been  compared  with  the  output  from 
the  backward-recurrence  method,  when  applied  to  the  same  potential.  Regardless  of  the 
potential's  shape,  it  has  been  shown  that  the  backward-recurrence  method  produces  the  same 
output  as  the  transfer  matrix  method,  when  the  potential  does  not  cause  the  transfer  matrix 
method  to  become  instable. 

The  backward-recurrence  method  shows  great  promise  in  the  numerical  solutions  of  T 
versus  E  curves  for  potential  energy  profiles  encountered  in  real  devices,  like  the  resonant 
tunneling  diode  potential  shown  in  Figure  7.  In  contrast,  the  transfer  matrix  method  proves  itself 
to  be  incapable  of  any  reasonable  accuracy  in  this  case.  Other  practical  potential  profiles  for 
nanoelectronic  devices  also  are  probably  within  the  capabilities  of  this  algorithm. 

The  backward-recurrence  method  for  calculating  transmission  coefficients  has  been  shown 
to  be  worth  further  development.  The  transfer  matrix  method,  while  limited  in  applicability,  can 
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act  as  a  worthwhile  benchmark,  against  which  future  versions  of  the  backward-recurrence  code 
may  be  tested. 
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APPENDIX  A. 
MATHEMATICAL  PROPERTIES  OF  TRANSFER  MATRICES 

A.         SYMMETRY  OF  THE  OVERALL  TRANSFER  MATRIX,  M 
1.    The  Example  of  One  Square  Barrier 

It  can  be  shown  that,  for  a  single  square  barrier  of  height  F0and  width  2a,  the 


overall  transfer  matrix,  M  - 


M(l,l)     M(l,2) 
M(2,l)     M(2,2) 


is  given  by 


M 


cosh  2Ka  h —  sinh  2Ka  \e2'ka 
2  ) 


ITJ 


sinh2Kor 


IT] 


t 


sinh2*ra 


IE 


cosh  2«or-  —  sinh  2kci  \e  2 

2 


where  k  =  J—^-(E  -  0)  is  the  wave  number  outside  the  barrier,  and  k  -  J—r-(V0  -  E)  is 
V  %  V  h" 

the  magnitude  of  the  purely  imaginary  wave  nurrtber  inside  the  barrier  {kbarner  =  ik  ). 
Also  in  this  expression  for  M  are  two  aggregate  constants,  rj  and  s ,  which  are  given  by 


77  = 


Kk     k 
(k     k 


yk         K 


Since,  in  this  construction,  k  and  k  are  purely  real,  s  and  77  are  real  as  well. 
(Merzbacher,  1961,  pp.  91-92) 

2.  The  General  Case 

The  symmetry  property  of  the  overall  transfer  matrix,  M,  is  obvious  from  its  form 
given  above.  M  has  the  form 
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M  = 


M(l,l)     M(l,2)' 
M(2,l)     M(2,2) 


a,  +b,i     a,  +bj 


•\  •  "\ 


a3+Z>3/'    a4+b4i 


where  ci2=ci3=0,  ai=a4,  bj=  -b4,  and  bi=  -b3.  In  other  words,  the  (1,2)  and  (2,1)  elements 
of  M  are  complex  conjugates  of  each  other,  and  are  both  purely  imaginary.  The  (1,1)  and 
(2,2)  elements  are  also  complex  conjugates  of  each  other,  but  each  has  both  real  and 
imaginary  components.  This  symmetry  is  not  limited  to  the  case  of  one  square  barrier; 
rather,  the  one  square  barrier  case  has  been  provided  as  an  illustrative  example. 


B.    DETERMINANT  OF  THE  OVERALL  TRANSFER  MATRIX,  M 

The  determinant  of  the  overall  transfer  matrix,  M,  is  identically  one.  This 
condition  is  true  for  all  potentials,  not  just  for  square  barriers.  This  property  can  be 
verified  in  the  M  given  for  the  square  barrier  example,  as  follows: 

V  ic-  \     ( 


M(l,l)     M(l,2) 
M(2,l)     M(2,2) 


cosh  2kci  h —  si  nh  2ko 

2 


cosh  2fca sinh  2kci 

2 


-l    T] 


sinh"  2fca 


V 


=  cosh2  2Ka  +  \ 

4       4 


sinh2  2kq  =  cosh2  2kq  -  sinh2  2kq  =  1. 


C.         SYMMETRY  AND  PROGRESSION  OF  INTERFACE  TRANSFER 
MATRICES,  N 

The  interface  transfer  matrices  have  different  symmetry  than  does  the  overall 
transfer  matrix,  M.  In  fact,  when  the  numerical  version  of  the  transfer  matrix  solution 
begins  to  degrade,  it  is  these  properties  of  the  interface  transfer  matrices  which  are 
involved.  The  following  is  output  from  the  transfer  matrix  program  for  a  condition 
known  to  cause  numerical  instability,  taken  from  Appendix  B,  section  B.3.  The  output  is 


for  a  sequence  of  square  barriers  whose  up-steps  occur  at  x=5, 1 5,25. . .  and  whose  down- 
steps  happen  at  x=  10,20,30... 

BARRIERS  =  10 


BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 
E=0. 1 96624  eV 

When  x=5. 000000,  N  is 
-1.3361 13e-01  +-8.990499e-02  j 
-1.3361 13e-01  +8.990499e-02  j 

When  x=l 0.000000,  N  is 
-1.799812e-01  +-1.479897e+01  j 
8.299747e-02  +8.161 762e-02  j  ' 

When  x- 15. 000000,  N  is 
-7.773521e-03  +-1.198140e-02  j 
-7.773521e-03  +1.198140e-02j 

When  x=20.000000,  N  is 
-6.730543e+01  +-1.527079e+02  j 
9.609614e-03  +3.772008e-03  j 

When  x=25. 000000,  N  is 

-2. 1 74860e-04  +- 1 .2478 1 5e-03  j 

-2.174860e-04+1.247815e-03j 

When  x=30.000000,  N  is 

-1. 37338  le+03  +-1.286365e+03  j 

9.151 546e-04  +-2.647852e-05  j 

When  x=3  5. 000000,  N  is 
2.565397e-05  +-1.093629e-04j 
2.565397e-05  +1.093629e-04  j 

When  x-40.000000,  N  is 
-1.993362e+04  +-7. 270083  e+03  j 
7.373917e-05  +-3.398770e-05  j 

When  x=45. 000000,  N  is 


■1.783441e+00  +  3.417415e-01j 
•1.783441e+00  +  -3.417415e-01  j 


-1.799812e-01  +  1.479897e+01  j 
8.299747e-02  +  -8.161762e-02j 


-2.001004e+01  + -4.341330e+00  j 
-2.001004e+01  +  4.341330e+00j 


-6.730543e+01  +  1.527079e+02j 
9.6096 14e-03  +  -3.772008e-03  j 


1.883595e+02  +  -1.3351 19e+02j 
1.883595e+02+  1.3351 19e+02j 


-1.373381e+03  +  1.286365e+03  j 
9.151 546e-04  +  2.647852e-05  j 


1.363410e+03  + -2.217759e+03  j 
•1.363410e+03  +  2.217759e+03  j 


-1.993362e+04  +  7.270083e+03  j 
7.373917e-05  +  3.398770e-05  j 


39 


5.896356e-06  +-8.029830e-06  j 
5.896356e-06  +8.029830e-06  j 

When  x=50.000000,  N  is 
-2.389105e+05  +1.273344e+04  j 
4.833746e-06  +-5.337304e-06  j 

When  x=5 5. 000000,  N  is 
7.602998e-07  +-4.500235e-07  j 
7.602998e-07  +4.500235e-07  j 

When  x=60.000000,  N  is 
-2.421817e+06  + 1.1 8850  le+06  j 
2.087238e-07  +-6.035369e-07  j 

When  x=65. 000000,  N  is 
7.767776e-08  +-1.027182e-08  j 
7.767776e-08  +1.027182e-08  j 

When  x=70.000000,  N  is 
-1.986518e+07  +2.303671e+07  j 
-3.961913e-09  +-5.649675e-08  j 

When  x=75. 000000,  N  is 
6.6943  3  5e-09  +1.863506e-09  j 
6.694335e-09  +-1.863506e-09  j 

When  x=80.000000,  N  is 
-1.041891e+08  +3.267910e+08  j 
-2.288082e-09  +-4.471323e-09  j 

When  x=85. 000000,  N  is 
4.813270e-10  +3.848460e-10  j 
4.813270e-10  +-3.848460e-10  j 

When  x=90.000000,  N  is 
3.643030e+08  +3.850378e+09  j 
-3.42171  le-10  +-2.852034e-10  j 

When  x=95. 000000,  N  is 
2.588327e-ll  +4.813621e-l  1  j 
2.588327e-l  1  +-4.813621e-l  1  j 

When  x=100.000000,  N  is 
2.080449e+10  +3.832754e+10  j 
-3.783386e-ll  +-1.136730e-ll  j 


-4.335545e+03  + -2.903269e+04  j 
-4.335545e+03  +  2.903269e+04  j 


-2.389105e+05  +  -1.273344e+04  j 
4.833746e-06  +  5.337304e-06  j 


8.340682e+04  + -3.203 155e+05j 
8.340682e+04  +  3.203 155e+05  j 


-2.421817e+06  +  -1.1 8850 le+06  j 
2.087238e-07  +  6.035369e-07  j 


2.281527e+06  +  -2.953693e+06  j 
2.281527e+06  +  2.953693e+06  j 


■1.986518e+07  + -2.30367  le+07j 
-3.961913e-09  +  5.649675e-08  j 


3.672609e+07  +  -2.054894e+07  j 
3.672609e+07  +  2.054894e+07  j 


-1.041891e+08  +  -3.267910e+08j 
-2.288082e-09  +  4.471323e-09  j 


4.718089e+08  + -5.07493  le+07j 
4.718089e+08  +  5.07493  le+07  j 


3.643030e+08  +  -3.850378e+09  j 
-3.42171  le-10  +  2. 852034e-10j 


5.118289e+09+  1.559857e+09j 
5.1 18289e+09  +  -1.559857e+09  j 


2.080449e+10  +  -3.832754e+10j 
-3.783386e-ll  +  1.136730e-l  1  j 
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Overall,  M  is 

-6.842685e+04  +5.333000e+05  j        4.766 148e+05  +  -2.488562e+05  j 
4.766148e+05  +2.488562e+05  j  -6.842685e+04  +  -5.333000e+05  j 

Broke  due  to  numerical  inaccuracy  @  E=0. 196624  eV 

The  program  output  shows  the  symmetry  and  progression  of  the  N  matrices.  At 

an  interface  where  a  step  increase  in  potential  occurs  (up-step),  the  form  of  N  is 


NuP 


c  +  di    a  +  bi 
c-di    a-  bi 


On  the  other  hand,  at  an  interface  where  a  step  decrease  in  potential  occurs  (down-step), 
N  has  the  form 


N 

down 


a  +  bi     a-  bi 

c  +  di     c  -  di 


where,  in  each  case,  a,  b,  c,  and  d  are  real  constants.  Note  that  the  Nup  and  Ndown  matrices 
have  distinctly  different  symmetry.  It  is  interesting  that  the  matrix  product  of  all  of  the  N 
matrices  has  diagonal  symmetry. 

Note  that,  as  x  increases,  the  magnitudes  of  a  and  b,  for  both  the  Nup  and  Ndown 
matrices,  get  very  large,  while  the  magnitudes  of  c  and  d  get  very  small.  At  the  point  that 
either  type  of  N  matrix'  c±di  elements  approach  the  numerical  limitations  of  double- 
precision  floating-point  numbers,  the  determinant  of  the  overall  transfer  matrix,  M, 
begins  to  diverge  from  one.  When  this  happens,  the  transfer  matrix  method  breaks  down, 
as  it  has  in  the  above  example. 
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APPENDIX  B. 

TRANSFER  MATRIX  METHOD  NUMERICAL  INSTABILITY 
AS  SEEN  IN  PROGRAM  OUTPUT 


A.         SINGLE-PRECISION  FLOATING-POINT  NUMBERS  USED 


Barrier  Height  Varied 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BAR  WIDTH  =  5.000000  nm 

Vo  =  0. 100000  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  5.000000  nm 

Vo  =  0.200000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 010702  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  5.000000  nm 

Vo  =  0.500000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 173308  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  5.000000  nm 

Vo=  1.000000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 680638  eV 


2.  Barrier  Width  Varied 


43 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 013768  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  6.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0.0361 12  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  7.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 072364  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  8.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E^0. 098584  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  =  9.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 1 18420  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BARWIDTH  -  10.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 161284  eV 
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3.  Number  of  Barriers  Varied 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  1 

BAR  WIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0.013768  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  2 

BARWIDTH  -  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 147604  eV 


%  ACC  withfloat 

%  withfloat 

BARRIERS  =  5 

BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0.221248  eV 


B.         DOUBLE-PRECISION  FLOATING-POINT  NUMBERS  USED 


1.  Barrier  Height  Varied 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  5.000000  nm 
Vo  =  0. 100000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  5.000000  nm 
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Vo  =  0.200000  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  1 

BAR  WIDTH  =  5.000000  nm 

Vo- 0.500000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  5.000000  nm 
Vo=  1.000000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  5.000000  nm 
Vo  =  2.000000  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  1 

BARWIDTH  =  5.000000  nm 

Vo  =  5.000000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 800680  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  1 

BARWIDTH  =  5.000000  nm 

Vo=  10.000000  eV 

Broke  due  to  numerical  inaccuracy  @  E=5. 629874  eV 


2.  Barrier  Width  Varied 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
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BAR  WIDTH  =  5.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  6.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  7.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  -  8.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  9.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  10.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  1 

BARWIDTH  =  20.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 003964  eV 
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%  ACC  withdouble 

%  withdouble 

BARRIERS  =  1 

BARWIDTH  =  30.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 107932  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  1 

BARWIDTH  =  40.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 161512  eV 

%  ACC  withdouble 

%  withdouble 

BARRIERS  -  1 

BARWIDTH  =  50.000000  nm 

Vo  -  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 187504  eV 


3.  Number  of  Barriers  Varied 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  1 
BARWIDTH  =  5.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  2 
BARWIDTH  =  5.000000  nm 
Vo  =  0.230000  eV 


%  ACC  withdouble 
%  withdouble 
BARRIERS  =  5 
BARWIDTH  -  5.000000  nm 
Vo  =  0.230000  eV 
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Broke  due  to  numerical  inaccuracy  @  E=0. 046372  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  10 

BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 196624  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  15 

BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0.2 15776  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  20 

BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 220792  eV 


%  ACC  withdouble 

%  withdouble 

BARRIERS  =  50 

BARWIDTH  =  5.000000  nm 

Vo  =  0.230000  eV 

Broke  due  to  numerical  inaccuracy  @  E=0. 227176  eV 
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APPENDIX  C. 
CODE  FOR  C  PROGRAMS  WRITTEN  FOR  THIS  THESIS 

A.         TRANSFER  MATRIX  METHOD  (methodl.c) 

/*  Francis  E.  Spencer  III  */ 
/*  Thesis.  Summer  1997*/ 
/*  Prof.  Luscombe  */ 

/*  Inclusions  */ 
#include  <stdio.h> 
#include  <math.h> 
#include  <stdlib.h> 

/*  Definitions  */ 

#define  BARRIERS  1  /*  dimensionless  */ 

#define  BARPTS  10  /*  dimensionless  */ 

#define  BAR  WIDTH  5.0  /*  nm  */ 

#define  MAXX  ((2.0*BARWIDTH*BARRIERS)+BARWIDTH)  /*  nm  */ 

#defme  XPTS  ((2*BARPTS*BARRIERS)+BARPTS)  /*  dimensionless  */ 

#define  EPTS  10000  /*  dimensionless  */ 

#define  EFFMASS  0.067  /*  dimensionless  */ 

#define  H20VER2M  0.0381  /*  eV-nm2  */ 

#define  C  (EFFMASS/H20VER2M)  /*  l/(eV-nm2)  */ 

#define  Vo  0.23  /*  eV  */ 

#define  ELEFT  0.000 1  /*  eV  */ 

#define  ERIGHT  (Vo-0.0001)         /*  eV  */ 

/*  Function  Prototypes  */ 

void  make V( void); 

struct  complex  addc(struct  complex  a  struct  complex  b); 

struct  complex  subc(struct  complex  a.  struct  complex  b); 

struct  complex  mulc(struct  complex  a,  struct  complex  b); 

struct  complex  divc(struct  complex  a  struct  complex  b); 

struct  complex  expc(struct  complex  a); 

double  absc(struct  complex  a); 

struct  complex  Nl  lcomplex(struct  complex  kminus,  struct  complex  kplus.  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N12complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N21complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N22complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  determinantcomplex(struct  complex  Mil.  struct  complex  Ml 2,  struct  complex  M21,  struct 

complex  M22); 

/*  Global  Variable  Definitions  */ 
struct  complex 

{ 

double  real; 
double  imag; 
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}; 

double  V[XPTS1; 
struct  complex  k[XPTS]; 

/*  Body  of  Program  follows: */ 

void  main(void) 

/*  This  function  controls  program  execution  */ 

{ 

int  ec.xc; 

double  E.epts.xpts.T; 

struct  complex 
newMl  l,newM12,newM21.newM22,oldMl  l,oldM12.oldM21,oldM22,Nl  l,N12,N21,N22,x,kratio,Mdet; 

x.imag=0.0;  /*  x  is  a  purely  real  number  */ 

epts=EPTS; 

xpts=XPTS; 

makeV(); 

for(ec=0.0;ec<epts;ec+=1.0) 

{ 

E=ERIGHT-(ERlGHT-ELEFT)*(ec/epts);    /*  DOWN -counting  through  E  values  */ 

oldM  1 1 .  imag=oldM  1 2.  imag=oldM2 1 .  imag=oldM22 .  imag=0 . 0; 

oldMll.real=oldM22.real=1.0; 

oldM12.real=oldM21.real=0.0;        /*  "old"  M  initially  an  identity  matrix  */ 

for(xc=0;xc<XPTS;xc+-f)  /*  initialize  k  */ 

{ 

if(E>=V[xc])        /*  k  is  purelv  real  */ 

{ 

k[xc].real=sqrt(C*(E-V|xc])); 

k[xc].imag=0.0; 
}/*endifE>=V*/ 
else  /*  k  is  purely  imaginary  */ 

{ 

k[xc].real=0.0; 

k[xc].imag=sqrt(C*(V[xc]-E)); 
}/*  end  else  E<V  */ 
}/*  end  for  xc  (init  k)  */ 

for(xc- 1  ;xc<XPTS;xc++) 

{ 

x.real=(xc/xpts)*MAXX; 

if((k[xcl.real!=k[xc-l].real)||(k[xc].imag!=k[xc-l].imag)) 

{ 

kratio=divc(k[xc],k[xc-l]); 

Nl  1=N1  lcomplex(k[xc-l],k[xc],kratio,x), 
N12=N12complex(k[xc-l],k[xc],kratio,x); 
N21=N21complex(k[xc-l].k[xcj.kratio.x); 
N22=N22complex(k[xc-l],k[xcl.kratio,x), 

newMll=addc(mulc(oldMll,Nll),mulc(oldM12,N21)); 
newM  1 2=addc(mulc(oldMl  1  ,N  12).mulc(oldM  1 2.N22)); 
newM2 1  =addc(mulc(oldM2 1  ,N  1 1  ).mulc(oldM22,N2 1 )); 
newM22=addc(mulc(oldM21.N12),mulc(oldM22.N22)); 

oldMll=newMll; 
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oldM12=newM12 
oldM21=newM21 
oldM22=newM22 
}/*endifk[.\c]  ...  */ 
}/*  end  for  xc  */ 

Mdet=detenninantcomplex(newMll.newM12,newM21,newM22); 

if((Mdet.real>1.0001)||(Mdet.real<0.9999)) 

{ 

printf("Broke  due  to  numerical  inaccuracy  @  E=%f  eV\n\n".E); 

break; 
}  /*  THE  "SAFETY"  IS  ON  */ 

else 
{ 


} 


T=sqrt(l,0/(absc(newMl  l)*absc(newMl  1))); 
pnntf("%.  12f\t%.  12f  \n'\E,T); 


}/*  end  for  ec  */ 
}/*  end  MAIN  */ 

void  makeV(void) 

/*  This  function  initializes  die  potential  energy  array,  V  */ 

{ 

int  xcount,y; 
for(xcount=0;xcount<XPTS;xcount++) 

{ 

y=xcount/BARPTS; 
if((y%2)==0) 

V[xcount]=0.0; 
else 

V[xcount]=Vo; 
}/*  end  for  xcount  */ 
}/*  end  MAKEV  */ 

struct  complex  addc(struct  complex  a,  struct  complex  b) 
/*  This  function  adds  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  sum; 

sum. real=a.real+b. real; 

sum.imag=a.imag+b.imag; 

return(sum); 
}/*  end  ADDC  */ 

struct  complex  subc( struct  complex  a.  struct  complex  b) 

/*  This  function  subtracts  complex  number  b  from  complex  number  a  */ 

{ 

struct  complex  difference; 

difference.real=a.real-b.real; 

difference .  imag=a.  imag-b .  imag; 

return(difference) ; 
}/*  end  SUBC  */ 

struct  complex  mulc(  struct  complex  a.  struct  complex  b) 

/*  This  function  multiplies  two  complex  numbers  passed  to  it  */ 
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{ 

struct  complex  product; 
product.  real=a.real*b.real-a.imag*b.imag; 
product.imag=a.real*b.imag+a.imag*b.real; 
return(product); 
}/*  end  MULC  */ 

struct  complex  divc(struct  complex  a.  struct  complex  b) 

/*  This  function  divides  complex  number  a  bv  complex  number  b  */ 

{ 

struct  complex  quotient; 

double  denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.  real=(a.real*b.real+a.imag*b.imag)/denom; 

quotient.  imag=(b.real*a.imag-a.real*b.imag)/denom; 

return(quotient); 
}/*  end  DIVC  */ 

struct  complex  expc(struct  complex  a) 

/*  This  function  computes  the  exponential  of  a  complex  number  */ 

{ 

struct  complex  exponential; 

exponential. real=exp(a.real)*cos(a.imag); 

exponential. imag=exp(a.real)*sin(a.imag); 

return(exponential); 
}/*  end  EXPC  */ 

double  absc(struct  complex  a) 

/*  This  function  returns  the  magnitude  of  complex  number  a  */ 

{ 

return(sqrt(a.real*a.real+a.imag*a.imag)); 
}/*  end  ABSC  */ 

struct  complex  Nl  lcomplex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  Nl  1  of  the  interface  */ 

{ 

struct  complex  one,  j,  N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(addc(one,kraUo),expc(mulc0.mulc(x,subc(kplus.kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*endNHCOMPLEX*/ 

struct  complex  N12complex(struct  complex  kminus.  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N 12  of  the  interface  */ 

{ 

struct  complex  one,  minusj,  N; 
one.real=1.0; 
one.imag=0.0; 
minusj.  real=0.0; 
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minusj.imag=-1.0; 

N=mulc(subc(one.kratio).expc(muJc(minusj,muJc(x,addc(kplus.kminus))))); 
N.real*=0.5; 
N.imag*=0.5; 
return(N); 
}/*endN12COMPLEX*/ 

struct  complex  N21complex(struct  complex  kminus.  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N2 1  of  the  interface  */ 

{ 

struct  complex  one,  j,  N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(subc(one,kratio),expc(mulc(j.mulc(x,addc(kplus.kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*  end  N21  COMPLEX*/ 

struct  complex  N22complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N22  of  the  interface  */ 

{ 

struct  complex  one,  minusj.  N; 

one.real=1.0; 

one.imag=0.0; 

minusj.  real=0.0; 

minusj. imag=-l  .0; 

N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 
}/*  end  N22COMPLEX  */ 

struct  complex  determinantcomplex(struct  complex  Mil,  struct  complex  Ml  2,  struct  complex  M21,  struct 

complex  M22) 

/*  This  function  calculates  the  determinant  of  the  transfer  matrix  as  a  check  for  accuracy  */ 

{ 

return(subc(mulc(Ml  l.M22),mulc(M2 1,M12))); 
}/*  end  DETERMINANTCOMPLEX  */ 


B.         BACKWARD-RECURRENCE  METHOD  (method2.c) 

/*  Francis  E.  Spencer  III  */ 
/*  Thesis,  Summer  1997*/ 
/*  Prof.  Luscombe  */ 

/*  Inclusions  */ 
^include  <stdio.h> 
#include  <math.h> 
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#include  <stdlib.h> 


/*  Definitions  */ 

#define  BARRIERS  1  /*  dimensionless  */ 

#define  BARPTS  5000  /*  dimensionless  */ 

/*  NOTE:  set  BARPTS  so  that  there  are  0.01  run  per  point  (BARWIDTH/BARPTS)=(  1/100)  */ 

#define  BARWIDTH  5.0  /*  nm  */ 

#define  MAXX  ((2.0*BARWIDTH*BARRIERS)+BAR  WIDTH)  /*  nm  */ 

#define  XPTS  ((2*BARPTS*BARRIERS)+BARPTS)  /*  dimensionless  */ 

#define  EPTS  10000  /*  dimensionless  */ 

#define  EFFMASS  0.067  /*  dimensionless  */ 

#define  H20VER2M  0.0381  /*  eV-nm2  */ 

#define  C  (EFFMASS/H20VER2M)  /*  l/(eV-nm2)  */ 

#define  Vo  0.23  .  /*  eV  */ 

#define  VL  0.0  /*  eV  */ 

#define  VM  0.0  /*  eV  */ 

#define  VR  0.0  /*  eV  */ 

#define  ELEFT  0.000 1  /*  eV  */ 

#defineERIGHT(Vo-0.0001)         /*  eV  */ 

/*  Function  Prototypes  */ 
void  makeV(void); 

struct  complex  addc(struct  complex  a,  struct  complex  b); 
struct  complex  subc(struct  complex  a,  struct  complex  b); 
struct  complex  mulc(struct  complex  a,  struct  complex  b); 
struct  complex  divc(struct  complex  a,  struct  complex  b); 
struct  complex  expc(struct  complex  a); 
double  absc(struct  complex  a); 

/*  Global  Variable  Definitions  */ 
struct  complex 

{ 

double  real; 
double  imag; 

}; 

double  V[XPTS]; 

/*  Body  of  Program  follows: */ 

void  main(void) 

/*  This  function  controls  program  execution  */ 

{ 

int  ec,xc; 

double  E,epts,xpts,barpts.T,delta.kl,kr; 

struct  complex  j.minusl.x,rn.rnplusl.ro,bn[XPTS],temp,temp2; 

x.imag=0.0;  /*  x  is  a  purely  real  number  */ 

epts-EPTS; 

xpts=XPTS; 

barpts=BARPTS; 

delta=BARWIDTH/barpts: 

j.real=0.0; 

j.imag=1.0; 

minus  l.real=- 1.0; 

minus  l.imag=0.0; 

makeV(); 
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for(ec=0.0;ec<epts;ec+=  1.0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts);    /*  DOWN-counting  througli  E  values  */ 
for(xc=0;xc<XPTS;xc++)  /*  initialize  bn*/ 

{ 

bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 

bn(xc].imag=0.0; 
}/*  end  for  xc  (init  k)  */ 

kl=sqrt(C*(E-VL)); 
kr=sqrt(C*(E-VR)); 
mplus  1  .real=cos(kr*delta); 
rnplus  1 .  imag=sin(kr*  delta); 


left  flat  zone  */ 


for(xc=(XPTS-BARPTS);xc>=BARPTS;xc«)  /*  edge  of  right  flat  zone  to  edge  of 

{ 

rn=divc(minusl,addc(bn[xc],rnplusl)); 

rnplus  1  real=rn.real; 

rnplus  1 .  imag=rn.  imag; 
}/*  end  for  xc  */ 

ro.real=rn.real; 
ro.imag=rn.imag; 

temp.real=0.0; 
temp.  imag=kl  *delta; 
temp2=mulc(ro,expc(temp)); 

T=sqrt((sui(kl*delta)/sui(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 
2.0*temp2.real+absc(ro)*absc(ro))); 

printf("%.  12f\t%.  12f\t%.  12f\t%.  12f\n",E,T); 

}/*  end  for  ec  */ 
}/*  end  MAIN  */ 

void  makeV(void) 

/*  This  function  initializes  the  potential  energy  array,  V  (square  barriers)  */ 

{ 

int  xcount,y; 

f or(xcount=0 ;  xcount<XPTS ;  xcount++) 

{ 

v=xcount/BARPTS; 
if((y%2)==0) 

{ 

if(xcount<BARPTS)  V[xcount]=VL;  I*  LEFT  flat  zone  */ 

else  if(xcount>(XPTS-BARPTS))  V[xcount]=VR;      /*  RIGHT  flat  zone  */ 
else  V[xcount]=VM;  /*  flat  zone  between  barriers  */ 

}/*  end  if  y  */ 
else 

V[xcount]=Vo;     /*  inside  a  barrier  */ 
}/*  end  for  xcount  */ 
}/*  end  MAKEV  */ 

struct  complex  addc(struct  complex  a,  struct  complex  b) 
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/*  This  function  adds  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  sum; 

sum.  real=a.real+b.  real; 

sum. imag=a.imag+b. imag; 

return(sum); 
}/*  end  ADDC  */ 

struct  complex  subc(struct  complex  a,  struct  complex  b) 

/*  This  function  subtracts  complex  number  b  from  complex  number  a  */ 

{ 

struct  complex  difference; 

difference.real=a.real-b.real; 

difference.  imag=a.  imag-b.  imag; 

return(difference) ; 
}/*  end  SUBC  */ 

struct  complex  mulc(struct  complex  a.  struct  complex  b) 

/*  This  function  multiplies  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  product; 

product.  real=a.real*b.real-a.imag*b.  imag; 

product.  imag=a.real*b.imag+a.imag*b.  real; 

return(product); 
}/*  end  MULC  */ 

struct  complex  divc(struct  complex  a,  struct  complex  b) 

/*  This  function  divides  complex  number  a  bv  complex  number  b  */ 

{ 

struct  complex  quotient; 

double  denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 

quotient.imag=(b.real*a.imag-a.real*b.imag)/denom; 

return(quotient); 
}/*  end  DIVC  */ 

struct  complex  expc(struct  complex  a) 

/*  This  function  computes  the  exponential  of  a  complex  number  */ 

{ 

struct  complex  exponential; 

exponential. real=exp(a.real)*cos(a.imag); 

exponential. imag=exp(a.real)*sin(a.  imag); 

return(exponential); 
\l*  end  EXPC  */ 

double  absc(struct  complex  a) 

/*  This  function  returns  the  magnitude  of  complex  number  a  */ 

retum(sqrt(a.real*a.real+a.imag*a.imag)); 
}/*  end  ABSC  */ 


C.         TRANSFER  MATRIX  METHOD  COMPARED  TO  ANALYTIC 
METHOD  (mlwithref.c) 
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/*  Francis  E.  Spencer  III  */ 
/*  Thesis.  Summer  1997*/ 
/*  Prof.  Luscombe  */ 

/*  Inclusions  */ 
#include  <stdio.h> 
#include  <math.h> 
#include  <stdlib.h> 

/*  Definitions  */ 

#define  BARRIERS  1  /*  dimensionless  */ 

#define  BARPTS  10  /*  dimensionless  */ 

#define  BAR  WIDTH  5.0  /*  nm  */ 

#define  MAXX  ((2.0*BARVvTDTH*BARRffiRS)+BARVVTDTH)  /*  nm  */ 

#define  XPTS  ((2*BARPTS*BARRIERS)+BARPTS)  /*  dimensionless  */ 

#define  EPTS  10000  /*  dimensionless  */ 

#define  EFFMASS  0.067  /*  dimensionless  */ 

#define  H20VER2M  0.038 1  /*  eV-nm2  */ 

#define  C  (EFFMASS/H20VER2M)  /*  l/(eV-nm2)  */ 

#define  Vo  0.23  /*  eV  */ 

#define  ELEFT  0.000 1  /*  eV  */ 

#define  ERIGHT  (5.0*Vo)  /*  eV  */ 

/*  Function  Prototypes  */ 

void  makeV(void); 

struct  complex  addc(struct  complex  a.  struct  complex  b); 

struct  complex  subc(struct  complex  a,  struct  complex  b); 

struct  complex  mulc(struct  complex  a,  struct  complex  b); 

struct  complex  divc(struct  complex  a,  struct  complex  b); 

struct  complex  expc(struct  complex  a); 

double  absc( struct  complex  a); 

struct  complex  Nl  lcomplex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N12complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N21complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N22complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  determinantcomplex(struct  complex  Mil,  struct  complex  Ml 2,  struct  complex  M21,  struct 

complex  M22); 

double  Tref  finder(void); 

/*  Global  Variable  Definitions  */ 
struct  complex 

{ 

double  real; 
double  imag; 

}; 

double  V[XPTS]; 
struct  complex  k[XPTS]; 

/*  Body  of  Program  follows: */ 

void  main(void) 

/*  This  function  controls  program  execution  */ 
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{ 

int  ecxc: 

double  E.epts,xpts,T; 
struct  complex 
newMll.newM12,newM21,newM22.oldMll,oldM12,oldM21,oldM22,Nll.N12,N21,N22,x,kratio,Mdet. 
x.imag=0.0;  /*  x  is  a  purely  real  number  */ 

epts=EPTS; 
xpts=XPTS; 
makeV(); 

for(ec=0.0;ec<epts;ec+=  1.0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts);    /*  DOWN-counting  through  E  values  */ 

oldM  1 1 .  imag=oldM  1 2.  imag=oldM2 1  .imag=oldM22.imag=0.0; 

oldMl  l.real=oldM22.real=1.0; 

oldM12.real=oldM21.real=0.0;        /*  "old"  M  initially  an  identity  matrix  */ 

for(xc=0;xc<XPTS;xc++)  /*  initialize  k  */ 

{ 

if(E>=V[xc])        /*  k  is  purely  real  */ 

{ 

k[xc].real=sqrt(C*(E-V[xc])); 

k[xc).imag=0.0; 
}/*endifE>=V*/ 
else  /*  k  is  purely  imaginary  */ 

{ 

k[xc].real=0.0; 

k[xcj.imag=sqrt(C*(V[xc]-E)); 
}/*  end  else  E<V  */ 
}/*  end  for  xc  (init  k)  */ 

for(xc=  1  ;xc<XPTS;xc++) 

{ 

x.real=(xc/xpts)*MAXX; 

if((k[xc).real!=k[xc-l].real)||(k[xc].imag!=k[xc-l].imag)) 

{ 

kratio=divc(k[xc],k[xc- 1  ]); 

Nll=Nllcomplex(k[xc-l],k[xc],kratio.x); 
N12=N12complex(k[xc-l],k[xc],kratio,x); 
N21=N21complex(k[xc-l],k[xc],kratio,x); 
N22=N22complex(k[xc-l],k[xc],kratio,x); 

newMll=addc(mulc(oldMll,Nll),mulc(oldM12,N21)); 
newM12=addc(mulc(oldMll,N12).mulc(oldM12,N22)); 
newM2 1  =addc(mulc(oldM2 1  ,N  1 1  ),mulc(oldM22,N2 1 )); 
newM22=addc(mulc(oldM21,N12).mulc(oldM22,N22)): 

oldMll=newMll; 
oldM12=newM12, 
oldM21=ne\vM21; 
oldM22=newM22; 
}/*endifk|xcl  ...  */ 
}/*  end  for  xc  */ 
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/*Mdet=deteniunantcomplex(newMll,newM12,newM21,newM22);*/ 

/*if((Mdet.real>1.0001)||(Mdet.real<0.9999)) 

{ 

printf("Broke  due  to  numerical  inaccuracy  @  E=%f  eV\n\n",E); 

break; 
}*/         /*  THE  "SAFETY"  IS  OFF  */  /*  restore  Mdet  line  above,  also,  to  turn  it  on  again  */ 
/*else 
{*/ 


/*}*/ 


T=sqrt(1.0/(absc(newMl  l)*absc(newMl  1))); 

printf("%.  12f\t%.  12f\t%.  12f\t%.  12An",E,T,Tref_finder(),(T-Tref_finder())); 


}/*  end  for  ec  */ 
}/*  end  MAIN  */ 

void  makeV(void) 

/*  This  function  initializes  the  potential  energy  array,  V  */ 

{ 

int  xcounty; 
for(xcount=0;xcount<XPTS;xcount-H-) 

{ 

y=xcount/BARPTS; 
if((y%2)== 0) 

V[xcount]=0.0; 
else 

V[xcount]=Vo; 
}/*  end  for  xcount  */ 
}/*  end  MAKEV  */ 

struct  complex  addc(struct  complex  a,  struct  complex  b) 
/*  This  function  adds  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  sum; 

sum.real=a.real+b.real; 

sum.imag=a.imag+b.  imag; 

return(sum); 
}/*  end  ADDC  */ 

struct  complex  subc(struct  complex  a,  struct  complex  b) 

/*  This  function  subtracts  complex  number  b  from  complex  number  a  */ 

{ 

struct  complex  difference; 

difference .  real=a.  real -b .  real ; 

difference.  imag=a.  imag-b.  imag; 

return(difference); 
}/*  end  SUBC  */ 

struct  complex  mulc(struct  complex  a,  struct  complex  b) 

/*  This  function  multiplies  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  product; 

product,  real  =ra.real*b.real-a.imag*b.  imag; 

product.imag=a.real*b.imag+a.imag*b.real; 

return(product); 
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}/*  end  MULC  */ 

struct  complex  divc(struct  complex  a,  struct  complex  b) 

/*  This  function  divides  complex  number  a  by  complex  number  b  */ 

{ 

struct  complex  quotient; 

double  denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 

quotient.  imag=(b.real*a.imag-a.real*b.imag)/denom, 

rerurn(quotient); 
}/*  end  DIVC  */ 

struct  complex  expc(struct  complex  a) 

/*  This  function  computes  the  exponential  of  a  complex  number  */ 

{ 

struct  complex  exponential; 

exponential.real=exp(a.real)*cos(a.imag); 

exponential. imag=exp(a.real)*sin(a.imag); 

return(exponential); 
}/*  end  EXPC  */ 

double  absc(struct  complex  a) 

/*  This  function  returns  the  magnitude  of  complex  number  a  */ 

{ 

return(sqrt(a.real*a.real+a.imag*a.imag)); 
}/*  end  ABSC  */ 

struct  complex  Nl  lcomplex( struct  complex  kminus.  struct  complex  kplus.  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N 1 1  of  the  interface  */ 

{ 

struct  complex  one,  j,  N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(addc(one.kratio),expc(mulc(j,mulc(x,subc(kplus.kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*  end  Nil  COMPLEX*/ 

struct  complex  N12complex( struct  complex  kminus.  struct  complex  kplus,  struct  complex  kratio.  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N 12  of  the  interface  */ 

{ 

struct  complex  one,  minusj,  N; 

one.real=1.0; 

one.imag=0.0; 

minusj.  real=0.0; 

minusj.  imag=- 1.0; 

N=mulc(subc(one,kratio),expc(mulc(ntinusj,mulc(x,addc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
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}/*  end  N12COMPLEX  */ 

struct  complex  N21complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio.  struct 

complex  x) 

/*  This  function  finds  die  matrix  element  N2 1  of  the  interface  */ 

{ 

struct  complex  one,  j,  N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=muIc(subc(one,kratio),expc(mulc(j,mulc(x,addc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*  end  N21  COMPLEX*/ 

struct  complex  N22complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N22  of  the  interface  */ 

{ 

struct  complex  one,  minusj,  N; 

one.real=1.0; 

one.imag=0.0; 

minusj.  real=0.0; 

minusj. imag=-l  .0; 

N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*  end  N22COMPLEX  */ 

struct  complex  determinantcomplex(struct  complex  Mil,  struct  complex  Ml  2,  struct  complex  M21,  struct 

complex  M22) 

/*  This  function  calculates  the  determinant  of  the  transfer  matrix  as  a  check  for  accuracy  */ 

{ 

return(subc(mulc(Mll,M22),mulc(M21,M12))); 
}/*  end  DETERMTNANTCOMPLEX  */ 

double  Treffinder(void) 

/*  This  function  computes  the  analytic  transmission  coefficient 

for  the  one-barrier  system  of  height  Vo  and  width  BAR  WIDTH  */ 

{ 

struct  complex  numdenom,kl,k2,two,four,barwidth,j,terml,term2; 

two.real=2.0; 

two.imag=0.0; 

four.real=4.0; 

four.imag=0.0; 

barwidth.real=BAR  WIDTH; 

barwidth.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

kl.real=k[0].real. 

kl.imag=k[0].imag; 

k2.real=k[BARPTS+l].real; 

k2.imag=k[B ARPTS+ 1  ]  imag; 
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num=mulc(mulc(mulc(expc(muJc(mulc(subc(k2,kl).barv\'idth),j)).k2),kl),four). 
terml=mulc(addc(kl.k2),addc(kl.k2)); 

term2=mulc(mulc(subc(k  1  .k2).subc(k  1  .k2)).expc(mulc(mulc(mulc(k2,banvidth)j).two))); 
denom=subc(tenn  1  .term2); 

return(sqrt(absc(divc(num.denom))*absc(divc(num.denom)))); 
}/*  end  TREFFINDER  */ 


D.         BACKWARD-RECURRENCE  METHOD  COMPARED  TO  ANALYTIC 
METHOD  (m2withref.c) 

/*  Francis  E.  Spencer  III  */ 
/*  Thesis.  Summer  1997*/ 
/*  Prof.  Luscombe  */ 

/*  Inclusions  */ 
#include  <stdio.h> 
#include  <math.h> 
#include  <stdlib.h> 

/*  Definitions  */ 

#define  BARRIERS  1  /*  dimensionless  */ 

#defme  BARPTS  500  /*  dimensionless  */ 

/*  NOTE:  set  BARPTS  so  that  there  are  0.01  run  per  point  (BAR WIDTH/BARPTS)=(  1/100)  */ 

#define  BAR  WIDTH  5.0  /*  run  */ 

#defme  MAXX  ((2.0*BARWIDTH*BARRIERS)+BAR  WIDTH)  /*  nm  */ 

#define  XPTS  ((2*BARPTS*BARRIERS)+BARPTS)  /*  dimensionless  */ 

#define  EPTS  10000  /*  dimensionless  */ 

#define  EFFMASS  0.067  /*  dimensionless  */ 

#define  H20VER2M  0.038 1  /*  eV-nm2  */ 

#define  C  (EFFMASS/H20VER2M)  /*  l/(eV-nm2)  */ 

#define  Vo  0.23  /*  eV  */ 

#defme  VL  0.0  /*  eV  */ 

#define  VM  0.0  /*  eV  */ 

#define  VR  0.0  I*  eV  */ 

#define  ELEFT  0.000 1  /*  e V  */ 

#define  ERIGHT  (5.0*Vo)  /*  eV  */ 

/*  Function  Prototypes  */ 

void  makeV(void); 

struct  complex  addc(struct  complex  a.  struct  complex  b); 

struct  complex  subc(struct  complex  a,  struct  complex  b); 

struct  complex  mulc( struct  complex  a.  struct  complex  b); 

struct  complex  divc(struct  complex  a.  struct  complex  b); 

struct  complex  expc(struct  complex  a); 

double  absc(struct  complex  a); 

double  Tref  finder(double  kl,  double  E); 

/*  Global  Vanable  Definitions  */ 
struct  complex 

{ 

double  real; 
double  imag; 

}; 

double  V[XPTS]; 
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/*  Body  of  Program  follows: */ 

void  main(void) 

/*  This  function  controls  program  execution  */ 

{ 

int  ec,xc; 

double  E.epts,xpts,barpts,T,delta.kl,kr; 

struct  complex  j, minus l.x.rrurnplusl,ro,bn[XPTS],temp,temp2; 

x.imag=0.0;  /*  x  is  a  purely  real  number  */ 

epts=EPTS; 

xpts-XPTS; 

barpts-BARPTS; 

delta=BARWIDTH/barpts; 

j.real=0.0; 

j.imag=1.0; 

minus  l.real=- 1.0; 

minus  l.imag=0.0; 

makeVO; 

for(ec=0.0;ec<epts;ec+=1.0) 

{ 

E-ERIGHT-(ERIGHT-ELEFT)*(ec/epts);    /*  DOWN-counting  through  E  values  */ 
for(xc=0;xc<XPTS;xc++)  /*  initialize  bn*/ 

{ 

bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 

bn[xc].imag=0.0; 
}/*  end  for  xc  (init  bn)  */ 

kl=sqrt(C*(E-VL)); 
kr=sqrt(C*(E-VR)); 
rnplus  1 .  real=cos(kr*  delta) ; 
rnplus  1 .  imag=sin(kr*delta); 

for(xc=(XPTS-BARPTS);xc>=B  ARPTS;xc--)  /*  edge  of  right  flat  zone  to  edge  of 

{ 


left  flat  zone  */ 


rn=divc(minusl,addc(bn[xc],rnplusl)); 
rnplus  1 .  real=rn.  real ; 
rnplus  1 .  imag=rn.  imag; 
}/*  end  for  xc  */ 

ro.real=rn.real; 
ro.imag=rn.imag; 

temp,  real =0.0; 
temp .  imag=kl  *  delta; 
temp2=mulc(ro,expc(temp)) ; 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 
2.0*temp2.real+absc(ro)*absc(ro))); 

printf("%.  12f\t%.  12f\t%.  12f\t%.  12f\n",E,T.Tref_finder(kl,E),(T-Tref_finder(kl,E))); 

}/*  end  for  ec  */ 
}/*  end  MAIN  */ 

void  makeV(void) 
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/*  This  function  initializes  the  potential  energy  array.  V  (square  barriers)  */ 

{ 

int  xcount,y; 
for(xcount=0.xcount<XPTS;xcount++) 

{ 

V=xcount/BARPTS; 
if((y%2)==0) 

{ 

if(xcount<BARPTS)  V[xcount]=VL;  /*  LEFT  flat  zone  */ 

else  if(xcount>(XPTS-BARPTS))  V[xcount]=VR;      /*  RIGHT  flat  zone  */ 
else  V[xcount]=VM;  /*  flat  zone  between  barriers  */ 

}/*  end  if  y*/ 
else 

V[xcount]=Vo;     /*  inside  a  barrier  */ 
}/*  end  for  xcount  */ 
}/*  end  MAKEV  */ 

struct  complex  addc(struct  complex  a,  struct  complex  b) 
/*  This  function  adds  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  sum; 

sum.real=a.real+b.real; 

sum.imag=a.imag+b.imag; 

return(sum); 
}/*  end  ADDC  */ 

struct  complex  subc(struct  complex  a,  struct  complex  b) 

/*  This  function  subtracts  complex  number  b  from  complex  number  a  */ 

{ 

struct  complex  difference; 

difference. real=a.real-b. real; 

difference.  imag=a.  imag-b.  imag; 

retum(  difference); 
}/*  end  SUBC  */ 

struct  complex  mulc(struct  complex  a,  struct  complex  b) 

/*  This  function  multiplies  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  product; 

product.  real=a.real*b.real-a.imag*b.  imag; 

product.  imag=a.  real  *b.  imag+a.imag*b.  real ; 

return(product); 
}/*  end  MULC  */ 

struct  complex  divc(struct  complex  a.  struct  complex  b) 

/*  This  function  divides  complex  number  a  by  complex  number  b  */ 

{ 

struct  complex  quotient; 

double  denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 

quotient.  imag=(b.real*a.imag-a.real*b.imag)/denom; 

return(quotient); 
}/*  end  DIVC  */ 

struct  complex  expc(struct  complex  a) 
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/*  This  function  computes  the  exponential  of  a  complex  number  */ 

{ 

struct  complex  exponential; 

exponential.  real=exp(a.real)*cos(a.imag); 

exponential .  imag=exp(a.  real )  *sin(a.  imag) ; 

retum(exponential); 
}/*  end  EXPC  */ 

double  absc(struct  complex  a) 

/*  This  function  returns  the  magnitude  of  complex  number  a  */ 

{ 

return(sqrt(a.real*a. real +a.imag*a. imag)); 
}/*  end  ABSC  */ 

double  Tref_finder(double  kl,  double  E) 

/*  This  function  computes  the  analytic  transmission  coefficient 

for  the  one-barrier  system  of  height  Vo  and  width  BAR  WIDTH  */ 

{ 

struct  complex  num,denom.kl,k2.two,four,barwidth,j,terml,term2; 

two.real=2.0; 

two.imag=0.0; 

four.real=4.0; 

four.imag=0.0; 

barwidth.real=BARWIDTH; 

barwidth .  imag=0 . 0 ; 

j.real=0.0; 

j.imag=1.0; 

kl.real=kl; 

kl.imag=0.0; 

if(E>= V[B ARPTS+ 1  ])      /*  k2  is  purely  real  */ 

{ 

k2.real=sqrt(C*(E-V[BARPTS+l])); 

k2.imag=0.0; 
}/*endifE>=V*/ 
else  /*  k2  is  purely  imaginary  */ 

{ 

k2.real=0.0; 

k2.imag=sqrt(C*(V[BARPTS+l]-E)); 
}/*  end  else  E<V  */ 
num=mulc(mulc(mulc(expc(mulc(mulc(subc(k2,kl),barwidth)j)),k2),kl),four); 
term  1  =mulc(addc(k  1  ,k2),addc(k  l,k2)); 

term2=mulc(mulc(subc(kl,k2),subc(kl,k2)),expc(mulc(mulc(mulc(k2,barwidth)j),two))); 
denom=subc(term  1  ,term2); 

return(sqrt(absc(divc(num,denom))*absc(divc(num,denom)))); 
}/*  end  TREF  FINDER  */ 


E.         TRANSFER  MATRIX  METHOD  APPLIED  TO  PARABOLIC  BARRIERS 

(mlpara.c) 

/*  Francis  E.  Spencer  III  */ 
/*  Thesis,  Summer  1997*/ 
/*  Prof.  Luscombe  */ 

/*  Inclusions  */ 
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#include  <stdio.h> 
#include  <math.h> 
#include  <stdlib.h> 

/*  Definitions  */ 

#define  BARRIERS  1  /*  dimensionless  */ 

#define  BARPTS  1000  /*  dimensionless  */ 

#define  BAR  WIDTH  5.0  /*  nm  */ 

#define  MAXX  ((2.0*BARWIDTH*BARRIERS)+BAR  WIDTH)  /*  nm  */ 

#define  XPTS  ((2*BARPTS*BARRIERS)+BARPTS)  /*  dimensionless  */ 

#define  EFTS  1000  /*  dimensionless  */ 

#define  EFFMASS  0.067  /*  dimensionless  */ 

#define  H20VER2M  0.0381  /*  eV-nm2  */ 

#define  C  (EFFMASS/H20VER2M)  I*  l/(eV-nm2)  */ 

#define  Vo  0.23  /*  eV  */ 

#define  ELEFT  0.001         .  I*  eV  */ 

#define  ERIGHT  (2.0*Vo)  /*  eV  */ 

/*  Function  Prototypes  */ 

void  makeV(void); 

struct  complex  addc(struct  complex  a  struct  complex  b); 

struct  complex  subc(struct  complex  a,  struct  complex  b); 

struct  complex  mulc(struct  complex  a.  struct  complex  b); 

struct  complex  divc( struct  complex  a,  struct  complex  b); 

struct  complex  expc(struct  complex  a). 

double  absc(struct  complex  a); 

struct  complex  Nl  lcomplex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N12complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N21complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x); 

struct  complex  N22complex(struct  complex  kminus,  struct  complex  kplus.  struct  complex  kratio,  struct 

complex  x); 

struct  complex  determinantcomplex(struct  complex  Mil.  struct  complex  Ml  2,  struct  complex  M21,  struct 

complex  M22); 

/*  Global  Variable  Definitions  */ 
struct  complex 

{ 

double  real; 
double  imag; 

}; 

double  V[XPTS]; 
struct  complex  k[XPTS]; 

/*  Body  of  Program  follows: */ 

void  main(void) 

/*  This  function  controls  program  execution  */ 

{ 

int  ec.xc; 

double  E.epts.xpts.T, 

struct  complex 
newMll,newM12.newM21.newM22.oldMll,oldM12.oldM21,oldM22.Nll,N12,N21,N22,x,kraUo.Mdet; 

x.imag=0.0;  /*  x  is  a  purely  real  number  */ 

epts-EPTS; 
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xpts=XPTS; 
makeVQ; 


for(ec=0;ec<EPTS;ec++) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); 

oldMl  1  imag=oldM12.imag=oldM2 1  .imag=oldM22.imag=0.0; 

oldMl  l.real=oldM22.real=1.0; 

oldM  12.real=oldM2 1  real=0.0; 

for(xc=0,xc<XPTS;xc-H-) 

{ 

if(E>=V[xc]) 

{ 


} 
else 

{ 


k[xc].real=sqrt(C*(E-V[xc])); 
k[xc].imag=0.0; 


k[xc].real=0.0; 
k[xc].imag=sqrt(C*(V[xc]-E)); 


} 
} 
for(xc=  1  ;xc<XPTS;xc++) 

{ 

x.real=(xc/xpts)*MAXX; 

if((k[xc].real!=k[xc-l].real)||(k[xc].imag!=k[xc-l].imag)) 

{ 

kratio=divc(k[xc],k[xc-l]); 

Nl  1=N1  lcomplex(k[xc-l],k[xc],kratio,x) 
N12=N12complex(k[xc-l],k[xc],kratio,x) 
N21=N21complex(k[xc-l],k[xcj,kratio,x) 
N22=N22complex(k[xc-l],k[xc],kratio,x) 

newM  1 1  =addc(mulc(oldM  1 1  ,N  1 1  ),mulc(oldM  1 2  ,N2 1 )) 
newM  1 2=addc(mulc(oldMl  1  ,N  12),mulc(oldM  1 2,N22)) 
newM2 1  =addc(muJc(oldM2 1  ,N  1 1  ),mulc(oldM22,N2 1 )) 
newM22-addc(mulc(oldM21,N12),mulc(oldM22,N22)) 

oldMl  l=newMll; 
oldM12=newM12; 
oldM21=newM21; 
oldM22=newM22; 


} 

Mdet=determinantcomplex(newM  1 1  .newM  1 2,newM2 1  ,newM22); 
if((Mdet.real>1.0001)||(Mdet.real<0.9999)) 

{ 

printf("Broke  due  to  numerical  inaccuracy  (a),  E=%f  eV\n\n",E); 
break; 

} 
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else 

{ 

T=sqrt(  1  0/(absc(newMl  1  )*absc(newM  1 1 ))); 
if(T>=0.001) 

pnntf("%.  12f\t%.  l2f\n",E,T); 

} 

} 
}/*  end  MAIN  */ 

void  makeV(void) 

/*  This  function  initializes  the  potential  energy  array,  V  */ 

{ 

int  xcount,xcount2,barcount', 

double  barpts,x.barVfBARPTS].tempV[XPTS],avg,pts; 

xcount=0, 

barpts-BARPTS; 

for(x=(-2.0*Vo);x<=(2.0*Vo);x+=(4.0*VoVbarpts) 

{ 

barV[xcount]=Vo-(x*x)/(4.0*Vo); 

xcount++:, 
}/*  end  for  x  */ 
for(xcount=0;xcount<BARPTS;xcount++) 

{ 

V[xcount)=0.0; 
}/*  end  for  xcount  */ 
for(barcount=0;barcount<BARRIERS;barcount-H-) 

{ 

xcount2=0, 
for(xcount-BARPTS+(barcount*2*BARPTS);xcount<(2*BARPTS)+(barcount*2*BARPTS);xco 
unt++) 

{ 

V[xcount]=barV[xcount2] ; 

xcount2++; 
}/*  end  for  xcount  */ 

for(xcount-(2*BARPTS)+(barcount*2*BARPTS);xcount<BARPTS+((barcount+l)*2*BARPTS); 
xcount++) 

{ 

V[xcount|=0.0; 
}/*  end  for  xcount  */ 
}/*  end  for  barcount  */ 

pts=BARPTS/ 100.0; 

for(xcount=0;xcount<XPTS;xcount++) 

{ 

avg=0.0; 
if((V(xcount]!=V[xcount-l])&&(xcount!=0)) 

{ 

for(xcount2=0;xcount2<(BARPTS/100);xcount2++) 

{ 

avg+=V[xcount+xcount2] ; 
}/*  end  for  xcount2  (create  avg)  */ 

avg=avg/pts. 


70 


for(xcount2=0;xcount2<(BARPTS/100);xcount2++) 
{ 

tempV[xcount+xcount2]=avg; 
}/*  end  for  xcount2  (write  to  tempV)  */ 

xcount+=((B  ARPTS/ 1 00)- 1 ); 
}/*  end  if  V  */ 

else        /*  "V+  =  V-"  */ 

tempV[xcount]=V[xcount]; 
}/*  end  for  xcount  (create  tempV)  */ 

for(xcount=0;xcount<XPTS;xcount++) 

{ 

V[xcount]=tempV[xcount]; 
}/*  end  for  xcount  (transfer  to  V[ xcount])  */ 
}/*  end  MAKEV  */ 

struct  complex  addc(struct  complex  a,  struct  complex  b) 
/*  This  function  adds  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  sum; 

sum. real=a.real+b. real; 

sum. imag=a.imag+b. imag; 

return(sum); 
}/*  end  ADDC  */ 

struct  complex  subc(struct  complex  a,  struct  complex  b) 

/*  This  function  subtracts  complex  number  b  from  complex  number  a  */ 

{ 

struct  complex  difference; 

difference .  real=a.  real  -b .  real ; 

difference.  imag=a.  imag-b.  imag; 

return(difference); 
}/*  end  SUBC  */ 

struct  complex  mulc(struct  complex  a.  struct  complex  b) 

/*  This  function  multiplies  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  product; 

product. real=a.real*b.real-a.imag*b. imag; 

product.  imag=a.real*b.imag+a.imag*b.  real; 

return(product); 
}/*  end  MULC  */ 

struct  complex  divc(struct  complex  a,  struct  complex  b) 

/*  This  function  divides  complex  number  a  by  complex  number  b  */ 

{ 

struct  complex  quotient; 

double  denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.  real=(a.real*b.real+a.imag*b.imag)/denom; 

quotient.irnag=(b.real*a.irnag-a.real*b.imag)/denom; 

return(quotient); 
}/*  end  DIVC  */ 
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struct  complex  expc(struct  complex  a) 

/*  This  function  computes  the  exponential  of  a  complex  number  */ 

{ 

struct  complex  exponential; 

exponential. reaJ=exp(a.real)*cos(a.imag); 

exponential. imag=exp(a.real)*sin(a.imag); 

retum(exponential); 
}/*  end  EXPC  */ 

double  absc(struct  complex  a) 

/*  This  function  returns  the  magnitude  of  complex  number  a  */ 

{ 

return(sqrt(a.real*a.real+a.imag*a.imag)); 
}/*  end  ABSC  */ 

struct  complex  Nl  lcomplex(struct  complex  kminus.  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  Nl  1  of  the  interface  */ 

{ 

struct  complex  one,  j,  N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(addc(one.kratio),expc(mulc(j,mulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*endNHCOMPLEX*/ 

struct  complex  N12complex(struct  complex  kminus,  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N 12  of  the  interface  */ 

{ 

struct  complex  one,  minusj,  N; 

one.real=1.0; 

one.imag=0.0; 

minusj.  real=0.0; 

minusj.  imag=-1.0, 

N=mulc(subc(one,kraUo),expc(mulc(minusj,mulc(x,addc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 
}/*endN12COMPLEX*/ 

struct  complex  N21complex( struct  complex  kminus.  struct  complex  kplus,  struct  complex  kratio.  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N2 1  of  the  interface  */ 

{ 

struct  complex  one,  j,  N; 

one.real=1.0; 

one.imag^O.O; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(subc(one.kratio),expc(mulc(j,mulc(x,addc(kplus,kminus))))); 
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N.real*=0.5; 
N.imag*=0.5; 
return(N); 
}/*  end N21  COMPLEX*/ 

struct  complex  N22complex( struct  complex  kminus.  struct  complex  kplus,  struct  complex  kratio,  struct 

complex  x) 

/*  This  function  finds  the  matrix  element  N22  of  the  interface  */ 

{ 

struct  complex  one.  minusj,  N; 

one.real=1.0; 

one.imag=0.0; 

minusj.  real^O.O; 

minusj.  imag^-l.O; 

N=mulc(addc(one,kratio).expc(mulc(minusj,mulc(x.subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

return(N); 
}/*  end  N22COMPLEX  */ 

struct  complex  determinantcomplex(struct  complex  Mil,  struct  complex  Ml 2,  struct  complex  M21,  struct 

complex  M22) 

I*  This  function  calculates  the  determinant  of  the  transfer  matrix  as  a  check  for  accuracy  */ 

{ 

return(subc(mulc(Mll,M22),mulc(M21.M12))); 
}/*  end  DETERMINANTCOMPLEX  */ 


F.         BACKWARD-RECURRENCE  METHOD  APPLIED  TO  RESONANT- 
TUNNELING  DIODE  (RTD)  POTENTIAL  (m2RTD.c) 

/*  Francis  E.  Spencer  III  */ 
/*  Thesis,  Summer  1997*/ 
/*  Prof.  Luscombe  */ 


/*  Inclusions  */ 

#include  <stdio.h> 

#include  <math.h> 

#include  <stdlib.h> 

/*  Definitions  */ 

#define  PI  3.141592654 

/*  dimensionless  */ 

#define  MAXX  50.0 

/*  nm  */ 

#define  XPTS  5000 

/*  dimensionless  */ 

#define  EPTS  1000 

/*  dimensionless  */ 

#define  EFFMASS  0.067 

/*  dimensionless  */ 

#define  H20VER2M  0.0381 

/*  eV-nm2  */ 

#define  C  (EFFMASS/H20VER2M)            /*  l/(eV-nm2)  */ 

#define  Vo  1.45 

/*  eV  */ 

#defme  VL  -0. 1 

/*  eV  */ 

#define  VR  -0.3 

/*  eV  */ 

#define  ELEFT  (VL+0.001) 

/*  eV  */ 

#define  ERIGHT  (Vo-0.001) 

/*  eV  */ 

/*  Function  Prototypes  */ 
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void  make V( void): 

struct  complex  addc(struct  complex  a.  struct  complex  b); 
struct  complex  subc(struct  complex  a.  struct  complex  b); 
struct  complex  mulc(struct  complex  a.  struct  complex  b); 
struct  complex  divc(struct  complex  a.  struct  complex  b); 
struct  complex  expc( struct  complex  a); 
double  absc(struct  complex  a): 
double  absval(double  x); 

/*  Global  Variable  Definitions  */ 
struct  complex 

{ 

double  real; 
double  imag; 

double  V[XPTS]; 

/*  Body  of  Program  follows: */ 

void  main(void) 

/*  This  function  controls  program  execution  */ 

{ 

int  ec.xc; 

double  E,epts,xpts,barpts,T,delta,kl,kr; 

struct  complex  j,minusl,x.rn.rnplusl.ro,bn[XPTS].temp,temp2; 

x.imag=0.0:         /*  x  is  a  purely  real  number  */ 

epts=EPTS; 

xpts=XPTS: 

delta=MAXX/xpts: 

j.real=0.0; 

j.imag=1.0; 

minusl.real=-1.0; 

minusl.imag=0.0; 

makeV(): 

for(ec=0.0;ec<epts;ec+=  1.0) 

{ 

E=ERlGHT-(ERIGHT-ELEFT)*(ec/epts):    /*  DOWN-counting  through  E  values  */ 
for(xc=0;xc<XPTS;xc++)  /*  initialize  bn*/ 

{ 

bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 

bn[xc].imag=0.0; 
}/*  end  for  xc  (init  bn)  */ 

kl=sqrt(C*(E-VL)); 

kr=sqrt(C*(E-VR)); 
rnplus  1  .real=cos(kr*delta); 
rnplus  1 .  imag=sin(kr*  delta); 

for(xc=(XPTS- 1  );xc>=0;xc~)         /*  edge  of  right  flat  zone  to  edge  of  left  flat  zone  */ 

{ 

rn=divc(  minus  l,addc(bn|xc|,  rnplus  1)); 

rnplus  l.real=rn.real; 

rnplus  1 .  imag=rn.  imag; 
}/*  end  for  xc  */ 

ro.real=rn.real; 
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ro.imag=rn.imag; 

temp.real=0.0; 
temp.imag=kl*delta; 
temp2=rnulc(ro,expc(teinp)); 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 
2.0*temp2.real+absc(ro)*absc(ro))). 

printf("%f\t%f\n",E,T); 

}/*  end  for  ec  */ 
}/*  end  MAIN  */ 

void  makeV(void) 

/*  This  function  initializes  the  potential  energy  array,  V  */ 

{ 

int  xcount,xcount2.PTS; 

double  barpts,x,xtwo.tempV[XPTS),avg,pts.xpts,last; 

xpts=XPTS; 

for(xcount=0;xcount<XPTS;xcount-H-) 

{ 

x=xcount*(PI/xpts)-PI/2.0; 

V[xcount]=(-0.  l*atan(x))-0.2; 
}/*  end  for  xcount  (create  the  AT  AN  shape)  */ 

for(xcount=0;xcount<XPTS;xcount++) 

{ 

x=xcount*(MAXX/xpts); 

xtwo=xcount*(Pl/xpts)-PI/2.0; 

if((x>=18.8)&&(x<19.0)) 

V[xcountl=((-0.1*atan(xtwo))-0.2)+(1.55/0.2)*(x-18.8); 
else  if((x>=19.0)«fe&(x<21.0)) 

V[xcount]=1.45-(0.05/2.0)*(x-19.0); 
else  if((x>=21.0)&«&(x<21.2)) 

V[xcount]=1.40-(1.55/0.2)*(x-21.0); 
else  if((x>=21.8)&&(x<22.0)) 

V[xcount]=((-0.1*atan(xtwo))-0.2)-(0.2/0.2)*(x-21.8); 
else  if((x>=22.0)&&(x<24.0)) 

V[xcount]=-0.35-(0.05/2.0)*(x-22.0); 
else  if((x>=24.0)&&(x<24.2)) 

V[xcount]=-0.40+(0.2/0.2)*(x-24.0); 
else  if((x>=24.8)&&(x<25.0)) 

V[xcountH(-0.1*atan(xtwo))-0.2)+(1.55/0.2)*(x-24.8); 
else  if((x>=25.0)«fe&(x<27.0)) 

V[xcount]=1.35-(0.05/2.0)*(x-25.0); 
else  if((x>=27.0)&&(x<27.2)) 

V[xcountl=1.30-(1.55/0.2)*(x-27.0); 
else;        /*  leave  V  unchanged  */ 

}/*  end  for  xcount  (add  the  3  peaks)  */ 
}/*  end  MAKEV  */ 

double  absval(double  x) 

/*  This  function  returns  the  absolute  value  of  a  double,  x.  */ 

{ 
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if(x<0)  x=-x; 
return(x); 
}/*  end  ABSVAL  */ 

struct  complex  addc(struct  complex  a.  struct  complex  b) 
/*  This  function  adds  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  sum; 

sum.real=a.real+b.real; 

sum.imag=a.imag+b.imag; 

return(sum); 
}/*  end  ADDC  */ 

struct  complex  subc(struct  complex  a,  struct  complex  b) 

/*  This  function  subtracts  complex  number  b  from  complex  number  a  */ 

{ 

struct  complex  difference; 

difference. real=a.real-b. real; 

difference.  imag=a.  imag-b.  imag; 

return(difference); 
}/*  end  SUBC  */ 

struct  complex  mulc(struct  complex  a,  struct  complex  b) 

/*  This  function  multiplies  two  complex  numbers  passed  to  it  */ 

{ 

struct  complex  product; 

product.real=a.real*b.real-a.imag*b.imag; 

product.imag=a.real*b.imag+a.imag*b.real; 

return(product); 
}/*  end  MULC  */ 

struct  complex  divc(struct  complex  a,  struct  complex  b) 

/*  This  function  divides  complex  number  a  by  complex  number  b  */ 

{ 

struct  complex  quotient; 

double  denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.  real=(a.  real  *b .  real+a.  imag*b.  imag)/denom; 

quotient.  imag=(b.real*a.imag-a.real*b.imag)/denom; 

return(quotient); 
}/*  end  DIVC  */ 

struct  complex  expc(struct  complex  a) 

/*  This  function  computes  the  exponential  of  a  complex  number  */ 

{ 

struct  complex  exponential; 

exponential. real=exp(a.real)*cos(a.imag); 

exponential. imag=exp(a.real)*sin(a.  imag); 

return(exponential); 
}/*  end  EXPC  */ 

double  absc(struct  complex  a) 

/*  This  function  returns  the  magnitude  of  complex  number  a  */ 

{ 

retum(sqrt(a.real*a. real+a. imag*a.imag));  }/*  end  ABSC  */ 
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APPENDIX  D. 
FIGURES 


1    Square  Barrier,    0.23    eV,    5,0   nm 
Method   1    (Transfer  Matrices) 
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Figure  1 .  Transfer  Matrix  Method  Applied  to  a  Single  Square  Barrier 
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1    Square   Barrier,    0.23    eV,    5.0   run 

Method   1    (Transfer  Matrices,     "Safety"    Off) 

minus  Analytic   Expression    (for   1    square  barrier) 
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Figure  2.  Difference  Between  Analytic  Solution  and  Transfer  Matrix  Solution,  for  the 

Case  of  One  Square  Barrier 
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Figure  3.  Potential  Energy  Profile  for  Three  Parabolic  Barriers 
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2    Parabolic    Barriers,       0.23    eV,     5.0    nm 
Method   1    (Transfer  Matrices) 
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Figure  4.  Transfer  Matrix  Method  Applied  to  a  Two-Parabolic-Barrier  Potential 
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3    Parabolic    Barriers,    0.23    eV,    5 . 0    nm 
Method   1     (Transfer  Matrices) 
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Figure  5.  Transfer  Matrix  Method  Applied  to  a  Three-Parabolic-Barrier  Potential 
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5  Parabolic  Barriers,  0.23  eV,  5 . 0  nm 
Method  1   (Transfer  Matrices) 
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Figure  6.  Transfer  Matrix  Method  Applied  to  a  Five-Parabolic-Barrier  Potential 
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Potential  Energy  versus  Distance  for  an 
AlAs/InGaAs/InAs  Resonant  Tunneling  Diode, 
as  in  Luscombe's  Article 
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Figure  7.  Potential  Energy  Profile  for  a  Resonant  Tunneling  Diode 
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RTD    Potential    Profile    (erroneous) 
TRANSFER    MATRIX    "SAFETY"     OFF 
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Figure  8.  Transfer  Matrix  Method  Applied  to  the  Resonant  Tunneling  Diode  Potential 
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1  Square  Barrier,  0.23  eV,  5.0  nm 
Method  2  (Backward  Recurrence) 
Compared  to  Analytic  Expression 
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Figure  9.  Backward-Recurrence  Method  Applied  to  the  Square  Barrier  of  Figure  1 
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1  Square  Barrier,  0.23  eV,  5.0  nm 
Method  2  (Backward  Recurrence) 
Minus  Analytic  Express xon 
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Figure  10.  Difference  Between  Analytic  Solution  and  Backward-Recurrence  Solution, 

for  the  Case  of  One  Square  Barrier 
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2    Parabolic    Barriers,       0.23    eV,    5.0    ran 
Method   2    (Backward   Recurrence) 
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Figure  1 1 .  Backward-Recurrence  Method  Applied  to  the  Two-Parabolic-Barrier 

Potential  of  Figure  4 
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3    Parabolic    Barriers,    0.23    eV,    5.0    nm 
Method   2     (Backward   Recurrence) 
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Figure  12.  Backward-Recurrence  Method  Applied  to  a  Three-Parabolic-Barrier  Potential 
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5  Parabolic  Barriers,  0.23  eV,  5.0  nm 
Method  2  (Backward  Recurrence) 
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Figure  13.  Backward-Recurrence  Method  Applied  to  a  Five-Parabolic-Barrier  Potential 
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5    Parabolic    Barriers,     0.23    eV,     5.0    run 
Method    2    (Backward   Recurrence) 
Using   10,000    Energy   Points 
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Figure  14.   Backward-Recurrence  Method  Applied  to  the  Five-Parabolic-Barrier  Potential 
of  Figure  13,  Using  Ten  Thousand  Energy  Points 
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RTD    Potential    Profile 
Using  Method    2    (backward   recurrence) 
(10, 000    E   points) 
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Figure  15.  Backward-Recurrence  Method  Applied  to  the  Resonant  Tunneling  Diode 

Potential 
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